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This  final  scientific  report  summarizes  progress  and  accon^lish- 
ments  under  contract  F44620-75-C-0009.  The  period  of  support  under  this 
contract  was  from  August  1,  1974  to  September  30,  1977.  'During  this  time 
the  investigation  focused  on  various  details  relating  to  the  nim)erical 
prediction  of  buoyant  vortex  wakes,  and  also  on  the  prediction  of  radiative 
decay  rates  of  atmospheric  waves.  A calculation  of  the  complete  vortex 
wake  including  buoyancy  was  not  conpleted  by  the  termination  date.  How- 
ever, fundamental  contributions  were  made  on  issues  relating  to  numerical 
false  viscosity  effects  at  large  Reynolds  number,  and  to  the  prediction  of 
the  descent  velocity  of  a closely  interacting  vortex  pair,  papers  were 
presented  on  these  topics.  These  contributions  are  described\in  the  pre- 
sent report.  A manuscript  for  publication  on  the  latter  study  is  now  in 
preparation.  The  work  on  atmospheric  radiative  decay  rate  was  the  subject 
of  an  oral  paper  and  has  also  been  published. 
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I.  General  Vortex  Wake  Considerations 

The  original  objective  of  the  present  investigation  concerned  the 
mechanism  by  which  large  supersonic  aircraft  flying  in  the  stably  strati- 
fied stratosphere  would  generate  atmospheric  internal  gravity  waves.  There 
is  no  question  that  such  a mechanism  exists  in  principle  from  two  sources: 
first,  the  vertical  motion  of  the  aircraft  vortex  wake;  second,  the  ultimate 
collapse  of  a well -mixed  wake  region  in  stratified  surroundings.  Both 

mechanisms  are  known  to  exist.  The  first  has  been  discussed  fran  two 

1 2 

essentially  different  points  of  view  by  Saffinan  and  Lissaman  et  al  , with 
contradictory  conclusions.  The  second  has  received  some  attention  due  to 
its  relevance  to  detection  of  submerged  objects  in  a stratified  oceein  based 
on  the  emitted  waves  from  wake  collapse^.  With  respect  to  aircraft,  prior 
work  did  not  permit  a ready  assessment  of  the  efficiency  of  generation  or 
the  magnitude,  persistence,  and  extent  of  the  wave  field  left  behind  in  the 
atmosphere.  All  that  can  be  concluded  with  certainty  from  prior  studies  is 
that  waves  will  certainly  be  generated.  If  their  anplitude  and  extent  is  or 
might  in  the  future  become  significant,  then  such  waves  would  be  of  concern 
both  in  atmospheric  diagnostics  dealing  with  presumably  naturally  occurring 
waves  and  possibly  also  with  respect  to  aircraft  operation. 

The  intent  of  this  study  was  to  avoid  the  need  for  the  arbitrary  and 
partly  contradictory  and  inconsistent  assumptions  which  characterize  previous 
theoretical  work  on  stratified  vortex  dynamics,  and  base  the  investigation 
on  numerical,  finite  difference  integrations  of  the  governing  partial 
differential  equations.  Such  an  approach  avoids  the  need  to  make  assump- 
tions regarding  vortex  wake  shape,  vorticity  entrainment  or  detrainment,  and 
vortex  generation  due  to  buoyancy,  but  rather  determines  these  phenomena  as 
part  of  the  solution. 
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Finite  difference  numerical  integrations  have  been  successfully 
applied  in  calculating  the  flow  field  of  rising  atmospheric  thermals, 
vdiich  are  analogous  to  this  study  in  involving  coupled  effects  of  vorticity 
and  buoyancy.  Similar  techniques  have  also  been  used  to  calculate  the 
generation  of  gravity  waves  from  the  collapse  of  a uniformly  mixed  region 
in  a stratified  atmosphere  ’ * . Therefore  there  is  every  reason  for  ex- 
pecting that  modifications  and  adaptations  of  known  techniques  will  lead 
to  numerical  predictions  fron  which  one  can  assess  the  magnitude  and 
propagation  characteristics  of  vortex  wake  generated  gravity  waves. 

At  the  present  termination  time  of  this  contract,  such  predictions 
have  not  been  achieved.  Among  the  various  features  of  the  calculation  re- 
quiring special  treatment,  two  were  particularly  troublesome  and  needed 
considerably  more  effort  than  had  originally  been  anticipated.  These  were 
deterioration  of  numerical  accuracy  for  a reasonable,  fixed  grid  size  with 
increasing  Reynolds  number,  and  means  of  keeping  the  computing  mesh  fixed 
on  the  vortex  wake  as  it  moves  vertically  at  a variable  and  unknown  velocity. 
Procedures  were  developed  to  deal  with  both  problems.  They  appear  to  be 
very  promising.  However,  the  opportunity  to  further  develop  and  apply  them 
to  the  problem  of  wave  generation  by  vortex  wakes  now  no  longer  exists. 

Both  the  problems  of  sufficient  accuracy  at  high  grid  Reynolds  numbers  and 
maintenance  of  a vorticity-fixed  coordinate  systan  are  quite  general  and 
not  only  of  interest  for  the  particular  problem  considered  here.  Because 
of  the  interest  and  importance  of  these  matters. the  approach  developed  under 
the  present  contract  for  dealing  with  them  is  described  in  this  report  in 
some  detail. 
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i A number  of  specific  vortex  dynamics  and  decay  predictions  were  1 

I f obtained,  both  before  and  after  incorporation  of  the  scheme  to  improve 

I accuracy  at  large  grid  Reynolds  number.  These  are  all  for  the  isothermal, 

I non-buoyant  case.  These  results  are  also  presented  in  this  report.  While  5 

* the  question  of  the  importance  of  a vortex  wake  as  a generator  of  gravity 

I 

waves  remains  unanswered,  it  still  appears  that  one  way  to  answer  it  is  to  | 

carry  to  conpletion  the  approach  taken  in  this  study.  I 

:l 

II.  Initial  Conditions  j 

During  an  earlier  period  of  AFOSR  sponsored  research,  a method  for  | 

the  numerical  finite  difference  integration  of  the  coupled  Navier-Stokes  j 

j 

(in  conservation  of  vorticity  form)  and  energy  equations  was  programmed  and 

7 

applied  to  a problem  involving  two-dimensional,  natural  convection  . This  j 

program  was  used  as  a base  frcmi  vdiich  the  calculation  me^od  for  the  present  i 

problem  evolved.  ; 

i 

Since  the  intent  was  to  calculate  the  dynamics  of  the  late,  far 

; 

downstream  stage  of  wake  development  and  not  the  early  formation  stage,  an 
initial  flow  and  temperature  field  representing  a buoyant  vortex  wake  is 

I 

needed  to  start  the  calculation.  This  should  be  representative  of  a vorticity  J 

field  which  is  distributed  over  the  wake  cross  section  and  not  unrealistically 
concentrated  into  small  vortex  "cores". 

The  choice  for  a starting  velocity  field  was  the  inviscid,  steady 

8 9 ■ 

state  vorticity  distribution  given  in  Lamb  and  also  discussed  by  Batchelor  . \ 

This  two-dimensional  flow  is  the  analog  of  the  well-known  Hill's  spherical 

vortex.  It  can  be  thought  of  as  representing  a steadily  downward  moving,  i 

f 

inviscid  vortex  wake  in  coordinates  moving  vertically  with  the  wake,  in  a I 

t 

i 

> 


L 


/: 

1: 
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plane  located  at  a large  fixed  distance  behind  an  aircraft. 

The  flow  is  given  by 

2VJt  (kr)sin9 

-Li — i , r<  a 


’P  = 


k jQ(ka) 


= - V(r  - ^ )sine 


r > a 


Here  is  a stream  function,  V is  a constant  reference  velocity,  r 
and  0 are  cylindrical  coordinates  with  0*0  pointing  vertically  downwards. 
Jq  and  are  Bessel  functions,  of  order  zero  and  one,  and  k is  a constant 
given  by  the  zero  of  which  defines  the  radius  of  a circle  which  confines 
the  region  of  distributed  vorticity.  Thus 


ka  * 3.83171 
JgCka)  = -0.402759 

The  stream  function  above  gives  an  inner  flow  whose  vorticity  r<  is 
proportional  to  itself,  through 

n * k^ip 

Therefore  the  vorticity  vanishes  on  the  bounding  circle  r * a,  while 
beyond  it  the  flow  is  irrotational . The  inner  flow  has  interior  stagnation 
points  separated  by  the  distance  0.961023a  where  vorticity  and  stream- 
function  have  a maximum.  Thus  it  represents  the  vorticity  distribution  of 
an  interacting  vortex  pair,  sinking  downwards  with  velocity  V,  in  a coordi- 
nate system  fixed  upon  this  vortex  pair.  This  velocity  is  given  by 

V = r/6. 8339a,  where  T is  the  circulation  about  half  the  flow,  i.e. 

00  00 

r = ndxdy.  There  are  no  infinite  velocities  anyvdiere.  Inner  and 

-00  0 

outer  flow  match  in  velocity,  vorticity,  and  shear.  The  inner  flow  is  an 
exact  solution  of  the  inviscid,  unsteady  equations.  The  action  of  viscosity 


i 


X 


m 


causes  the  interior  vorticity  to  decay  and  to  spread  beyond  the  initial 
confining  boundary  given  by  r = a. 

These  features  are  illustrated  in  Fig.  1.  The  upper  part  of  the 
figure  shows  streamlines  on  the  left,  contours  of  constant  vorticity  on 
the  right.  These  results  ccane  fran  our  computer  program;  they  represent 
the  theoretical  solution  just  described  but  obtained  from  a finite  difference 
integration  of  several  time  steps  after  an  initial  input  of  the  theoretical 
solution.  In  such  a starting  procedure,  viscosity  is  taken  as  zero  in  the 
calculation.  The  lower  part  of  the  figure  shows  the  horizontal  distribution 
of  vertical  velocity  v on  the  left,  and  the  vorticity  on  the  right.  These 
quantities  have  been  normalized  by 

V V 

Here  v'  and  Q'  are  dimensional,  and  v is  kinematic  viscosity.  H is 
the  size  of  the  square  conputing  region  comprising  half  the  flow  field,  and 
this  is  four  times  the  initial  radius  of  the  vortex  bubble.  The  solid  lines 
reoresent  the  initial,  inviscid  distributions. 

The  dotted  curves  indicate  the  result  of  a viscous  calculation 
carried  out  to  a later  time.  T.ie  expected  re-distribution  of  vorticity  and 

attenuation  of  vertical  velocity  is  evident.  A dimensionless  time  is  defined 

, . I vt 

by  t = . 

For  the  non- isothermal  case,  no  ccmiputations  were  made  but  an  analo- 
gous temperature  field  to  start  the  calculation  has  been  formulated. 

The  initial  velocity  field  satisfies  the  appropriate  steady  inviscid 


equation  of  motion 


^“53r‘  w “5>r" 


through 


= -k\  . 
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The  steady,  non-conducting  version  of  the  energy  equation  admits 
solutions  of  the  same  kind  through  the  more  general  form 
4)  = 

where  (J)  = + AT  + yy  is  the  sum  of  ambient  temperature  T^,  an  unknown 

tffTiperature  perturbation  AT  to  be  solved  for,  and  compressibility  in  the 
ambient  atmosphere  represented  by  the  adiabatic  lapse  rate  y.  A particular 
case  would  be  4)  * in  vdiich  lines  of  eqioal  potential  temperature 

would  be  identical  to  lines  of  constant  vorticity.  This  is  not  an  appropri 
ate  starting  solution,  since  it  corresponds  to  a vortex  pair,  one  of  which 
is  hotter  and  the  other  colder  than  ambient.  A more  appropriate  solution 
appears  to  be 

(j>  = . 

This  function  represents  hot  spots  centered  on  the  vortices  vanishing 
smoothly  at  the  initial  wake  boundary  r = a where  4;  = 0.  The  temperatur'j 
excess  scales  with  (j)Q.  This  function  would  serve  as  an  initial  input  to 
the  energy  equation,  with  4)q  as  a disposable  parameter. 

III.  Boundary  Conditions 

The  flow  illustrated  in  Fig.  1 has  the  property  that  far  from  the 
origin  the  velocity  asymptotically  becomes  a spatially  constant  vertical 
updraft  with  magnitude  V.  This  boundary  condition  cannot  be  enforced  on 
the  finite  boundaries  of  the  computing  region  which  surrounds  the  vortex 
pair,  since  an  upper  bound  exists  on  the  size  of  this  control  surface  due 
to  computing  cost  limitations.  Typically,  the  ccmiputing  region  cannot  be 
larger  than  an  order  of  magnitude  greater  than  a.  The  problem  is  quite 
common  and  occurs  vdienever  boundary  conditions  at  infinity  need  to  be 
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imposed  with  finite  difference  calculations  and  finite  budgets.  The 
following  scheme  was  developed  to  deal  with  this  problem. 

The  interior  vorticity  distribution,  whatever  it  may  be,  determines 
the  velocity  at  arbitrary  points  x^,  y^^  through  the  appropriate  Green's 
function: 


H X 


Ujj(t)  = - I JJ JJCx',y',t) 


0 0 


[CV^')^+Cyb-y')^]  Cyb-y')^] 


i 


i 


H X 


y^it)  = vet)  + i j^n(x',y',t) 


0 0 


[x^-x ' ^ - (yjj-y ' ) ^]  X ' dx 'dy ' 

[ (Xjj-x ' ) (yj^-y ' ) [xjj+x ' ) (yj^-y ' ) 2] 


Here  x'  and  y'  are  coordinates  of  arbitrary  interior  points  where 
the  vorticity  is  n.  Integration  over  all  such  points  contained  in  a 
rectangular  control  surface  specifies  the  velocity  at  some  point  on  this 
surface  (x^,y|j),  vdiere  vertical  and  horizontal  boundaries  are  given  by  H 
and  X.  The  above  expressions  have  been  incorporated  into  our  numerical 
program,  and  they  evaluate  boundary  velocities  at  any  time  based  on  the 
vorticity  distribution  which  has  been  determined  at  that  time . These  boundary 
velocities  are  then  related  to  boundary  stream- functions  \diich  are  used  to 
find  interior  velocities  used  to  go  on  in  time. 

This  Green's  function  formulation  for  velocity  boundary  conditions 
has  worked  very  well.  Accuracy  as  measured  by  comparing  the  numerically 
calculated  inviscid  flow  including  boundary  values  after  a few  time  steps 
with  the  exact  analytical  results  for  this  flow  is  satisfactory.  This 
formulation,  being  kinematic,  is  equally  valid  for  isothermal  and  buoyant 


1 

I 

i 


calculations . 
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In  a stratified  ambient  environment,  the  vortex  wake  will  appear 
to  move  through  a time-dependent  temperature  field  if  observed  in  a 
coordinate  systan  fixed  on  the  wake.  If  Tg  is  the  temperature  at  any 
height  far  from  the  wake,  then  at  a point  there  fixed  relative  to  the  wake 
one  will  see  a temporal  temperature  change  given  by 

3T 

sT- 

Here  y is  the  adiabatic  lapse  rate.  This  will  be  the  appropriate  tempera- 
ture boundary  condition.  This  time -varying  ten^erature  can  be  applied, 
without  serious  error,  on  the  boundary  of  the  ccmputing  region,  in  contrast 
to  the  case  for  velocities.  This  boundary  condition  introduces  yet  another 
role  for  V(t).  For  positive  stratification + y > oj,a  downward  moving 
vortex  wake  (V<0)  will  sense  an  ever  cooler  environment  and  become  ever  more 
buoyant.  Thus  V attains  great  dynamical  significance  for  the  non- isothermal 
case. 


rv.  Conventional  Finite  Difference  Method  Results 

! 

f 

The  finite  difference  scheme  initially  utilized  is  standard  and  only 
the  treatment  of  the  non-linear  convection  terms  requires  conment.  A 
particular  three-point,  non-central  operator  proposed  by  Torrance^^’^^  was 
I used.  It  employs  forward  or  backward  differences  depending  on  the  direction 

i of  the  local  velocity.  The  operator  has  the  property  that  it  satisfies  con- 

f 

servation  requirements  and  imposes  no  grid  size  restriction  to  prevent  nu- 
merical instability.  These  strong  advantages  need  to  be  balanced  against 

t 

the  disadvantage  that  the  resulting  accuracy  is  only  first  order  in  grid 
size,  in  contrast  to  second  order  accurate  three  point  central  differences. 
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One  way  to  view  this  matter  is  that  with  non -central  differences,  results 
can  be  obtained  vdiich  may  become  inaccurate  at  too  coarse  a grid;  with 
central  differences,  results  cannot  then  be  obtained  at  all  due  to  numeri- 
cal instabilities.  The  general  impression  regarding  the  resulting  finite 
difference  equations,  gathered  at  the  time  this  study  began  from  the 
literature  and  personal  communications,  was  that  while  formally  some  in- 
accuracy problems  were  to  be  expected  at  large  grid  Reynolds  number;  in 
practice  these  were  not  serious. 

Results  obtained  with  this  method  for  the  isothermal  wake  will  now 
be  briefly  described. 

Two  runs  were  made . The  values  of  the  Reynolds  number  were  5 and 

250.  This  Reynolds  number  is  based  on  the  initial  downward  vortex  wake 

velocity  and  the  initial  spacing  between  the  vortex  centers  as  given  by  the 

wing  span.  If  the  latter  is  50m,  and  the  downward  velocity  is  1 msec'^ 

(corresponding  to  an  aircraft  speed  of  300  rasec'^  with  typical  cruise  lift 

conditions) , then  the  corresponding  values  for  viscosity  are  10^  and 
3 2 -1 

2 X 10  cm  sec  , respectively.  The  first  represents  a large  value  typical 
of  large  scale  eddy  mixing  in  the  troposphere,  the  second  is  more  typical 
of  turbulence  generated  within  the  wake  itself. 

The  initial  flow  configuration  has  already  been  shown  in  Fig.  1. 

At  the  small  Reynolds  number,  streamlines  and  vorticity  100  time  steps  into 
the  calculation  are  shown  in  Fig.  2.  This  time  corresponds  to  a downward 
motion  of  the  wake  of  10%  of  its  original  diameter.  It  has  grown  sub- 
stantially larger,  but  without  strong  distortion  of  the  bounding  streamline. 
For  the  large  Reynolds  number,  240  time  steps  result  in  the  flow  shown  in 
Fig.  3.  By  this  time  the  wake  has  moved  downwards  just  under  1.5  initial 
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diameters.  It  is  now  significantly  distorted,  as  shown  both  by  the 
vorticity  pattern  and  bounding  streamline.  A coiqxjsite  illustration  show- 
ing the  progressive  distortion  of  the  bounding  streamline  as  a function  of 
time,  for  both  Reynolds  numbers,  is  given  by  Fig.  4,  where  the  dots  near 
the  center  indicate  the  outward  movement  of  the  vortex  "center".  This 
distance  is  used  as  a measure  of  wake  size.  The  wake  boundary  is  super- 
inqxjsed  on  vorticity  contours  at  the  largest  time  for  each  run,  in  Fig.  5. 
This  illustrates  the  beginning  development  of  a vortex  wake  when  the  Reynolds 
number  is  large. 

As  the  wake  grows  and  distorts,  its  downward  velocity  adso  changes. 
For  these  exanples  it  decreases,  so  that  the  wake  slows  down.  For  both 
csises,  the  runs  were  of  such  duration  that  the  wake  drift  velocity  de- 
creased to  about  50%  of  its  initial  value. 

A particular  representation  of  this  phenomenon  is  shown  in  Fig.  6. 

The  ordinate  in  both  plots  is  a special  combination  of  downward  velocity, 

12 

wake  size,  and  wake  in^julse  suggested  by  an  analysis  due  to  Nfeixworthy  . 

The  quantity  I represents  a measure  of  the  impulse  needed  to  generate  the 
motion  from  rest.  Here,  it  has  been  evaluated  from  I * (j)^»  with  b the 
instantaneous  vortex  spacing.  The  combination  of  variables  making  up  the 
ordinate  of  Figs.  8 and  9 is  one  which  is  predicted  in  reference  12  to  be 
linearly  proportional  to  time.  It  was  found  possible  to  determine  an 
effective  origin  in  our  examples  such  that,  very  nearly,  the  predicted 
linear  variation  is  observed. 

Two  aspects  of  Fig.  6 need  further  discussion.  First,  the  Max- 
12 

worthy  prediction  is  not  only  a linear  variation  with  time  of  the  ordinate 
of  Fig.  6,  but,  in  our  present  non-dimensional  variables,  also  a slope  vdiich 
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is  independent  of  Reynolds  number.  This  is  not  found  to  be  the  case  with 
I the  initial  numerical  calculations.  It  was  not  clear  whether  this  Reynolds 

number  effect  on  slope  is  just  a consequence  of  wake  distortion,  which 
would  influence  wake  drag  and  decay  and  is  not  included  in  Maxworthy 's 
f analysis,  or  numerical  errors. 

The  other  matter  concerns  the  obvious  starting  transients  visible 
in  the  data.  A modification  has  been  made  subsequent  to  the  generation  of 
the  data  shown,  in  which  a new  initial  step  ranoves  the  discrepancy  between 
initial  analytical  data  and  a slightly  inaccurate  version  of  this  same  data, 
which  is  produced  by  numerical  relaxation. 

V.  Modification  to  Reduce  False  Transport  Effects 

As  the  investigation  proceeded,  several  indications  began  to  suggest 
that  all  was  not  well  with  the  high  Reynolds  number  calculation.  First,  for 
this  condition  it  was  found  that  the  downward  drift  velocity  would  still  de- 
cay significantly  even  with  deletion  of  viscosity.  Second,  the  decay  was 
significantly  reduced  by  changing  the  computing  mesh  from  60x60  to  120x120. 
While  the  disagreement  regarding  the  Maxworthy  prediction  of  Reynolds  number 
independence  was  also  disquieting,  this  disagreement  was  somewhat  ambiguous 
because  the  latter  prediction  is  approximate  and  not  intended  for  Reynolds 
numbers  as  low  as  those  considered  here.  Finally,  other  assessments  regard- 
ing the  Torrance  convective  operator  appeared^^  which  gave  definite  indi- 
cation that  the  Reynolds  ntimber  of  250  calculation  was  probably  seriously 
contaminated  by  false  transport  effects.  For  the  coarse  mesh  (60x60),  our 
highest  grid  Reynolds  number  was  approximately  10,  which  should  now,  in 
light  of  present  insight,  be  considered  unacceptably  large.  It  was  realized 
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that  if  false  transport  effects  were  present,  these  wov'^'i  seriously  affect 
subsequent  gravity  wave  predictions,  due  to  their  expected  relatively  fine 
spatial  structure. 

The  following  scheme  was  designed  in  an  attempt  to  eliminate  first 
order  grid  size  errors  without  going  to  central  differences,  thereby  hope- 
fully avoiding  the  need  for  ever  decreasing  mesh  size  with  increasing 
Reynolds  number  to  maintain  numerical  stability.  To  my  knowledge,  this 
very  sinqjle  idea  has  not  been  explored  before. 

From  a Taylor  series  expansion  of  the  function  whose  derivative  is 
to  be  approximated,  one  can  obtain  for  the  point  x*x^,  where  u*u^: 


3US2 


3ufl 


Ax 


r 3u  3nl 


OCAx^) 


u.  1+U.  u.+u.  , 

VO.,  0 } >-!■  > 0 

2 ’ 2 


3ua  3un  ix  r ^ 3u  3!i1 ...  n-,,2 


U.^,+U.  u.+u.  1 

0,  ^ < 0 


Here  u is  the  horizontal  velocity  in  the  x direction,  and  H is  the 
vorticity.  The  subscript  FD  denotes  finite  difference,  DIF  denotes  differ- 
ential. These  particular  forms  result  from  the  special  difference  operators 
used  here^^’^^.  If  opposite  signs  occur  for  the  mean  velocities  at  adjacent 
gridpoints,  no  such  simple  expression  appears  to  be  derivable.  For  that  case, 
the  approximation  is  suggested  that  the  first  order  terms  should  be  neglected 
entirely.  Corresponding  expressions  can  be  derived  for  the  vertical  con- 
vective transport  terms. 
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The  equations  above  are  not  new.  The  first  order  difference  be- 
tween finite  difference  and  differential  derivatives  is  usually  called 
false  viscosity  when  u or  v is  constant,  here  they  are  sin5)ly  referred  to 
as  first  order  errors.  Let  all  such  terms  be  luin)ed  together  as  EAx. 

If  the  ranaining  terms  in  the  original  unsteady  differential 
equation  for  vorticity  are  similarly  expanded,  one  finds  that  the  differ- 
ence equation  which  is  to  be  solved  numerically,  which  has  the  form 

FD  equation  » 0 (1) 

becomes 

2 

DIF  equation  + EAx  + ^ -^  + O(Ax^)  = 0 (2) 

^ 3t^ 

The  interpretation  of  (2)  is  that,  to  0(Ax),  the  differential 
equation  is  not  satisfied  if  one  has  programmed  (1) . But  this  is  only  a 
question  of  book-keeping.  The  term  EAx  can  be  attached  to  (1)  rather  than 
(2)  if  one  considers  the  difference  equation  to  be 

FD  equation  * + EAx  (3) 

One  would  then  find 

2 

DIF  equation  + ^ §-4  + OCAx^)  = 0 (4) 

2 at^ 

The  second  order  error  in  time  is  found  to  be  negligible  for  the 
problem  dealt  with  here. 

Equation  (3)  can  be  interpreted  as  the  result  of  constructing  a 
finite  difference  approximation  not  of  the  real  differential  equation  but 
of  an  artificial  one  in  such  a way  that  the  finite  difference  version  of 
the  artificial  equation  is,  to  second  order  accuracy,  the  proper  differen- 
tial equation. 
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Since  the  quantity  EAx  involves  various  space  derivatives,  these 
are  known  or  can  be  operated  with  in  exactly  the  same  way  as  the 
differential  operators  already  occurring  in  (1)  or  the  left  side  of  (3) . 

A scheme  to  improve  accuracy  therefore  is  to  program  and  apply  (3)  rather 
than  (1) . This  modification  to  our  original  program  has  been  made  and 
tested  in  a limited  way,  with  the  following  results. 

For  an  inviscid,  unmodified  run  (all  viscous  terms  removed  from 

the  governing  equation) , one  finds  with  the  normal  60  x 60  grid  a flow 

decay  of  0.394%.  With  modification  it  is  0.005%.  In  this  example  the 

correct  decay  is  known;  it  is  no  decay  at  all.  The  modified  program 

does  better  by  a factor  of  79,  which  is  of  the  order  expected  in  going  to 
2 

Ax  accuracy  con5)ared  to  Ax. 

Some  tests  were  also  made  for  the  viscous,  large  Reynolds  number 
case  described  in  the  previous  section.  Here  comparisons  were  made  be- 
tween modified  and  unmodified  program,  both  with  a normal  grid  and  a grid 
half  that  size.  The  modification  produced  a decay  which  was  not  signifi- 
cantly affected  by  halving  the  grid  size.  This  decay  is  slower  than  that 
with  the  unmodified  program  and  the  original  grid  size,  and  differs  even 
more  from  the  lower  Reynolds  number  case  than  was  the  case  with  the  un- 
modified program. 

Examination  of  the  flow  field  predicted  at  the  large  Reynolds  number 
with  the  modified  program  shows  the  development  of  small  oscillations 
in  the  vorticity  contours.  This  is  illustrated  in  Fig.  7.  This 
configuration  exists  after  60  time  steps.  It  is  not  known  at  this  time 
if  this  result  is  real,  and  the  possibility  needs  to  be  examined  that  a 
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flow  instability  is  developing  at  this  fairly  large  Reynolds  number.  The 
most  obvious  test  is  to  repeat  the  calculations  with  a finer  mesh;  such 
calculations  have  not  yet  been  made.  No  such  oscillations  develop  with 
the  unmodified  program  and  enhanced  false  viscosity  effects,  or  at  the 
low  Reynolds  number. 


VI.  Theoretical  Analysis  of  Vertical  Wake  Velocity 

The  instantaneous  overall  vertical  velocity  V plays  several  roles. 
It  is  an  a priori  unknown  variable  in  time  and  one  of  the  main  character- 
istics of  the  vortex  wake  to  be  found.  It  is  also  needed  to  keep  the  com- 
puting grid  centered  on  the  vortex  system.  Without  such  an  arranganent 
the  region  of  interest  would  soon  escape  any  coiqiuting  grid  of  reasonable 
size  and  resolution.  Finally,  it  is  required  for  the  temperature  boundary 
condition,  as  described  in  Section  III. 

A nimerical  control  system  was  constructed  which  monitors  the  in- 
terior stagnation  point  of  the  flow  field  for  vertical  drift  of  the 
calculated  flow  with  respect  to  the  coordinate  system  supposedly  moving 
vith  the  vortex  system.  A correction  was  then  calculated  from  this  data 
and  applied  to  V for  the  succeeding  time  step  to  prevent  this  drift.  The 
correction  was  designed  to  so  alter  V,  at  each  time  step,  that  net  drift 

is  prevented,  to  a tolerance,  for  all  time.  The  resulting  V(t)  is  then 
the  instantaneous  vertical  downward  velocity  of  the  vortex  wake 
which  is  sought. 

Some  problems  were  encountered  in  predicting  the  evolution  of  V(t). 
An  initial,  crude  version  of  a control  system  for  drift  prevention  was  too 


insensitive  and  produced  no  control;  a second  version  was  much  too  sensi- 


tive  and  produced  unstable  oscillations.  Several  succeeding  formulations 
all  failed  to  danp  these  oscillations. 

Finally,  however,  a schane  was  found  which  seemed  to  be  both  stable 
and  effective  in  preventing  drift.  This  consisted  of  three  separate 
contributions  to  a corrective  velocity  which  are  proportional  to  stagna- 
tion point  displacement,  displacement  rate,  and  the  time  derivative  of 
the  latter.  Judicious  choice  for  the  coefficients  of  these  corrections  was 
able  to  produce  a non-oscillatory  V(t)  with  the  unmodified,  original  ccm- 
puter  program,  and  this  scheme  was  used  to  produce  the  results  described 
in  Section  IV.  However,  the  situation  still  was  not  very  satisfactory.  A 
rational  and  reliable  means  for  determining  the  coefficients  was  never 
found.  Further,  there  were  cases  in  which  oscillations  would  begin  again 
after  a considerable  computing  time  with  smooth  decay,  and  then  new  co- 
efficients would  have  to  be  found,  by  trial  and  error. 

When  the  conqjuter  program  included  the  modifications  described  in 
Section  V to  remove  false  transport  effects,  the  problem  of  oscillations 
in  the  predicted  V(t)  returned.  At  a Reynolds  number  of  250,  coefficients 
in  the  correction  to  V were  never  found  to  truly  stabilize  the  motion  of 
the  coordinate  system.  This  is  indicated  in  Fig.  8. 

In  the  isothermal  case  this  situation  is  unfortunate  but  not  one  to 
entirely  prevent  proper  interpretation  of  calculated  results.  In  effect, 
one  views  the  flow  field  as  one  might  a television  screen  picture  with 
oscillations  in  the  vertical  hold  control.  At  any  instance,  the  spatial 
picture  is  correct,  and  the  oscillation  is  an  annoyance.  Since  the  wake 
velocity  V should  almost  certainly  decay  in  a non-oscillatory  manner,  one 
can  even  still  extract  V from  some  averaging  process.  However,  with  a 
buoyant  wake  in  a stratified  atmosphere  one  of  the  primary  questions  is 
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vrfiether  the  entire  wake  will  oscillate,  and  how  nuch.  Then  the  situation 
just  described  is  entirely  unacceptable.  Therefore  a different  method 
for  determining  the  quantity  V(t)  had  to  be  found.  This  section  briefly 
describes  the  results  of  work  done  on  this  problem. 

The  approach  taken  was  to  search  for  a theoretical  basis  for 
determining  the  overall  velocity  of  a region  of  distributed  vorticity, 
allowing  for  viscous  and  non-steady  effects.  No  derivation  of  a theo- 
retical expression  for  velocity  was  found  in  the  literature,  nor  was  our 
attempt  at  this  successful. 
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Even  for  the  incon5)ressible , isothermal  case,  only  Saffinan  appears 
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to  have  seriously  addressed  the  problem.  Although  Saffinan 's  study  is 
applied  to  a viscous  vortex  ring  of  small  cross-section  relative  to 
diameter,  he  also  gives  a general  definition  of  vortex  system  velocity. 

This  velocity  is  defined  rather  than  derived,  but  he  gives  supporting  argu- 
ments for  the  definition.  Therefore  Saffiman's  definition,  as  well  as  some 
others  based  on  intuition  or  found  in  the  literature,  was  calibrated  on 
two  cases  for  suitability  of  incorporating  such  theoretical  expressions  into 
the  corqjuter  program.  This  calibration  leads  to  the  conclusion  that  Saff- 
man's  definition  appears  to  be  correct  and  suitable,  as  will  now  be  described. 

Consider  a general  vorticity  distribution  n(x,y)  vdiich  is  anti- 
symmetric about  the  vertical  y axis.  For  the  sp>ecial  case  of  our  starting 
solution,  there  is  also  symmetry  about  the  horizontal  x axis,  but  with  time 
this  disappears  under  the  action  of  viscosity.  In  this  analysis,  the 


variables  have  not  been  normalized. 


The  quantity  TJdxdy  is  preserved  for  all  time  in  any  two-dimen- 
sional boundary- free  viscous  flow;  unfortunately  here  this  quantity  is  zero. 
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It  can  be  shown  that  the  non-vanishing 
the  vertical  impulse  is  also  preserved.  This 


integral  corresponding  to 
is  given  by 


•00 


A similar  invariant  integral  corresponds  to  the  horizontal  impulse 
I^;  for  the  present  case  this  is  zero. 

Saffinan^^  argues  that  the  vertical  velocity  of  the  entire  vorticity 
distribution  is  given  by 


V » 
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fMxdy 


(6) 


For  the  presait  case,  this  reduces  to 

00 

V » J^xy  Ifdxdy  (7} 

^ -00 

The  time  derivative  represents  the  variation  with  time  of  the 

vorticity  distribution  in  a coordinate  system  fixed  in  space  (not  on  the 

vortex  system),  at  a fixed  position. 

The  dynamic  equation  governing  the  evolution  of  vorticity  can  now 
go 

be  used  to  evaluate  ^ . 

For  the  isothermal,  incompressible  case  (7)  can  be  written  as 

00 

v=  - JJ  37  -vV^jxydxdy  (8) 

It  can  be  shown,  after  some  manipulations,  that  for  the  present 
case  (8)  can  eventually  be  written  as 
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V = 


J'J yu)fMxciy 

-CO 

00 

JJ ^fidxdy 


C9) 


Various  terms  arising  from  integration  by  parts,  and  evaluated  far 
away  from  the  region  near  n,  vanish  because  of  the  anti -symmetrical  nature 
of  the  vorticity  distribution. 

Once  the  definition  of  vortex  system  velocity  given  by  Saffinan^^ 
has  been  manipulated  into  the  form  (9) , direct  comparison  with  other 
definitions  becomes  possible. 

Not  only  does  the  quantity 
This  suggests  the  definition 


“ 00 

J'J I^dxdy  vanish,  but  so  does  Jj  vi^dxdy. 


V = 0 


00  00 

J'J'  vf2  dxdy 


If 


(10) 
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Q dxdy 

-00  0 

Lo  and  Ting’^*'  claim  the  equivalent  of  (10)  to  be  the  velocity  of  the  vortex 
system.  A similar  form  is  used  by  Bilanin  et  al^^. 

The  fact  that  the  denominator  of  (10)  will  generally  be  time  de- 
pendent may  be  a worry  to  some,  and  an  obvious  remedy  is  to  define 


(11) 


This  intuitive  definition  is  recognized  as  the  first  term  of  (9) . 
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The  velocities  from  (9),  (10),  and  (11)  will  be  denoted  by  Vg, 
Vj^g,  and  It  is  straightforwaid  to  evaluate  them  for  the  steady, 

inviscid  distributed  vorticity  which  has  been  the  starting  flow  of  our 
numerical  computations.  It  should  be  remembered  that  all  integrands  are 
referred  to  coordinates  fixed  in  space. 

If  the  known  steady  velocity  in  this  case  is  Vg  (see  Section  II), 


then  one  finds 


'^9  - ^10  ' ^'^11  ’ ''O- 


In  other  words,  both  terms  of  (9)  are  necessary,  and  both  Saffman's 
and  the  more  conventional  definition  (10)  are  correct  for  this  inviscid  steady 


flow. 


A second  comparison  of  the  various  expressions  for  V is  with  the 


viscous,  unsteady  flow  configuration  of  two  counter-rotating  Lamb  vortices. 
The  evolution  of  this  configuration  was  also  calculated  by  Bilanin  et  al^^, 
using  a numerical  finite  difference  model  including  turbulence.  Here  the 
approach  is  laminar,  and  analytical  predictions  can  be  obtained. 

Let  the  vorticity  be 


\ = Irfvt 


exp  - 


r ^ 2 

(x-j)  +y 


Thus  two  vortices  of  opposite  sign  are  located  at  x = tj,  y=0, 
separated  in  x by  the  distance  b.  The  initial  circulation  about  each  is  ±r. 
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Whatever  the  separation,  early  in  time  the  vortices  do  not  interact  since 
their  viscous  cores  (av^)  are  much  smaller  than  b.  At  large  times  they 
vdll  strongly  interact,  each  attenuating  the  vorticity  of  the  other.  In 
evaluating  the  various  integrals,  at  each  time  t the  spatially  fixed 
coordinate  system  is  taken  with  the  x axis  coincident  with  the  instan- 
taneous vertical  vortex  center  position. 

2 

As  t Q,  b /4vt  -►  <»,  the  vortex  cores  become  6 functions  and  Vg  * 
v^j=  r/2Trb,  the  latter  quantity  being  the  correct,  known  steady 
velocity  for  that  case. 

2 

As  t -*■<»,  b /4vt  -*•  0,  one  finds 


1 rb 

S’  4irvt 


'^10-  -I’®! 


^ll“  I ^9 


(12) 

(13) 


In  both  the  present  examples  the  two  contributions  to  the  numerator 
of  (9J  were  found  to  be  equal,  explaining  vdiy  = y Vg.  This  is  probably 
due  to  the  special  nature  of  the  test  flows:  both  maintain  vertical  symmetry 
about  the  x axis,  without  the  trailing  tail  of  detrained  vorticity  which 
developed  in  our  numerical  COTqjutations  and  also  known  to  exist  in 
traveling  viscous  vortex  configurations  generally. 

The  finding  that  vortex  system  velocity  is  proportional  to  the  quan- 
tity rb/vt  is  believed  to  be  new.  It  is  characteristic  of  the  very  last 

stage  of  vortex  interaction,  and  differs  (but  only  slightly)  from  that  pre- 
12 

dieted  by  Maxworthy  for  an  earlier  stage  of  wake  decay.  To  facilitate 


jM 
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comparison,  the  present  prediction  in  the  dimensionless  variables  of 
Fig.  6 is 


-1 

V I a t. 


The  constant  of  proportionality  given  by  (12)  differs  from  that  in 
(13) . The  correct  theoretical  value  is  not  known.  All  one  can  say  is  that 
there  are  supporting  arguments  in  favor  of  Vg  while  none  seem  to  exist  for 


^10' 


Eqiiation  (9 J is  attx  .ctive  for  other  reasons  besides  a veneer  of 
theoretical  rigor.  The  integrands  in  both  numerator  and  denominator  tend 
to  zero  as  the  y axis  is  approached  (fJ,  x,  and  u vanish  there)  more  rapidly 
than  those  in  (10).  As  mentioned  before,  there  will  in  general  be  a vertical 
trail  of  detrained  vorticity  located  downstream  of  the  vorticity  distri- 
bution and  centered  about  the  y axis.  This  must  eventually  cross  one  of 
the  horizontal  boundaries  of  a vortex-fixed  ccmputing  region  and  be  lost  to 
the  computation,  thereby  introducing  errors.  Use  of  (9)  will  minimize  these. 

The  conclusion  of  this  part  of  the  investigation  is  that  for  the 
isothermal  case  equation  (9),  expressed  in  vortex  fixed  coordinates,  should 
be  incorporated  into  the  computer  program  and  used  as  a basis  for  adjustments 
of  the  coordinate  system  velocity.  For  the  general  variable  density  case. 


Saffinan’s  definition^^  of  velocity  might  still  serve  as  a base  for  deriving 


a generalization  to  (10) . This  expectation  is  based  on  the  fact  that  I 


still  retains  its  incompressible  form  in  the  Boussinesq  approximation,  as 
1 


shown  by  Saffinan 


1 
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VII.  Recanmendations  for  Future  Work 

1.  Equation  C9)  should  be  incorporated  into  the  computer  program 
as  a means  of  adjusting  the  velocity  of  the  vortex-fixed  coordinate 
system. 

2.  The  incompressible  form  of  C9)  should  be  generalized  to  include 
variable  density,  utilizing  the  Boussinesq  approximation. 

3.  The  scheme  for  eliminating  first  order  (in  grid  size)  computing 
errors  should  be  further  tested.  This  aspect  of  the  investigation  is 
potentially  the  most  valuable  and  useful,  as  it  extends  beyond  the  con- 
fines of  the  particular  problem  studied  here.  Tests  to  separate  real  from 
computational  instability  can  be  made  vdth  the  present  problem  through  a 
sequence  of  runs  at  progressively  smaller  grid  size,  and  also  by  constructing 
numerical  solutions  for  other  loiown,  stable  flows. 

4.  The  effectiveness  by  which  a strong  vortex  wake  in  a stable  at- 
mosphere generates  internal  gravity  waves  should  also  be  examined  from  a 
physical,  less  ''brute  force"  point  of  view, to  conqjlement  any  future  studies 
such  as  these.  The  present  study  has  not  led,  during  the  contracting  period, 
to  an  evaluation  of  this  mechanism.  At  this  time  it  is  neither  possible  to 
dismiss  it,  nor  to  advocate  it  as  a strong  potential  influence  on  the  upper 
atmosphere  or  on  aircraft  flying  there. 


27 


ACKNCWLEDGEMENTS 

The  progranming  and  numerical  work  carried  out  under  this  con- 
tract would  have  been  impossible  without  the  outstanding  abilities  of 
Ms.  S.  Yamaimra  during  the  early  phase  of  the  work  and,  subsequently, 
those  of  Mr.  Carl  Hutton.  I wish  to  acknowledge  the  courage,  support, 
and  understanding  of  the  contract  monitor,  Milton  Rogers. 


28 


References 

1.  Saffinann,  P.  G.,  "The  Motion  of  a Vortex  Pair  in  a Stratified 
Atmosphere",  Studies  in  Applied  Mathematics  2,  107-119,  1972. 

2.  Lissaman,  P.B.S.,  Crow,  S.C.,  McCready,  P.B.  Jr.,  Tombach,  I.H. , 
and  Bate,  E.R.  Jr. ,"  Aircraft  Vortex  Wake  Descent  and  Decay  under  Real 
Atmospheric  Effects",  AeroVironment  Report  No.  FAA-RD-73-120,  AD  771-311/8, 
1973. 

3.  Schooley,  A.  H.,  and  Bighes,  B.A. , "An  Experimental  and  Theo- 
retical Study  of  Internal  Waves  Generated  by  the  Collapse  of  a Two- 
Dimensional  Mixed  Region  in  a Density  Gradient",  J.  Fluid  Mech.  1, 

159-175,  1972. 

4.  Wessel,  W.R. , "Numerical  Study  of  the  Collapse  of  a Perturbation 
in  an  Infinite  Density  Stratified  Fluid",  Phys.  Fluids,  1^,  1211,  171-176, 
1969. 

5.  Young,  J.  A.,  and  Hirt,  C.  W. , "Numerical  Calculations  of  Internal 

Wave  Motions",  J.  Fluid  Mech.,  2,  256-276,  1972. 

6.  Orlanski,  I.,  and  Ross,  B.,  "Numerical  Simulation  of  the  Generation 
and  Breaking  of  Internal  Gravity  Waves",  J.G.R. , 7^,  36,  8808-8826,  1973. 

7.  Traugott,  S.  C.,  and  Yamamura,  S.H. , "Buoyant  Convection  in 

Radiating  Horizontal  Fluid  Layers",  AIAA  J.,  1,  88-94,  1975. 

8.  Lamb,  H.,  "Hydrodynamics",  6th  Edition,  Dover  Publication,  p.  245, 

1932. 

9.  Batchelor,  G.  K. , "An  Introduction  to  Fluid  Dynamics",  Cambridge 
University  Press,  p.  535,  1967. 

10.  Torrance,  K.E.,  "Ccmiparison  of  Finite  Difference  Computations  of 
Natural  Convection",  J.  Res.  N.B.S.,  72B,  4,  281-300,  1968. 


29 


L 

t 

11.  Torrance,  K.  E.,  and  Rockett,  J.  A. , "Numerical  Study  of 

f Natural  Convection",  J.  Fluid  Mech. , 36,  1,  33-54,  1969. 

> — 

12.  Maxworthy,  T.,  "The  Structure  and  Stability  of  Vortex  Rings", 

J.  Fluid  Mech.,  1,  15-32,  1972. 

^ 13.  DeVahl  Davis,  G.,  and  Mallinson,  C.  D.,  "An  Evaluation  of 

Upwind  and  Central  Difference  Approximations  by  a Study  of  Recirculating 
Flow",  Canputers  and  Fluids,  29-43,  1976. 

14.  Saffinan,  P.  G.,  "The  Velocity  of  Viscous  Vortex  Rings",  Studies 
in  Appl.  Math.,  XLIX,  4,  371-380,  1970. 

15.  Lo,  K.  C.,  and  Ting,  L.,  "Studies  of  the  Merging  of  Vortices", 
Phys.  Fluids,  19,  6,  912-913,  1976. 


16.  Bilanin,  A.  J.,  Teske,  M.  E. , and  Williamson,  C.  G.,  "Vortex 
Interactions  and  Decay  in  Aircraft  Wakes",  AIAA  J. , 2,  250-260,  1977. 


Fig.  6 Vertical  Wake  Velocity  Decay. 


81 


IME  STEPS 


w 


t 


APPENDIX:  PROGRAM  LISTING 
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1 


t 


I 

f 

1 


//HAMCJ126  JOB  ( 3 1 66.  OR  A3 1 6600  . CO  .MML 1 3 . 25 1 7 1 1 1 *60  ) t ••  VORTEX  PAIR", 

//  CLASS*G, 

//  MSGLevEL=(l,n ,CONO»(7,UTI 

//  EXEC  FORTXCUG,  

//  XREF=NOXREF, 

//  G0TIME=5.'  ^ 

//  GOCORE  = 150k  _ 

//FORT.SYSIN  00  * ~ ■ "■  

C 

C PROGRAM  TO  TEST  AOOl TiON" OP  VlSCSsTfV  TO  INVisCiO  VORTEX  PAIR 

COMMQN/BOX/DT.TIHE.OX.DY.H.PSI (33,63) , VOHT ( 33. 63 ,2 > . U ( 33 ,63) .V (33. 

1631  .UMAX. VMAX.CAPU.OELUU)  ,K0(vT  ! 

COMMON/INOEX/IX.IY  i 

COMMON/BOX l/  TcPxn3, 63)  .FCFY  (33,631 

C ,VC1T(33.63) ,VC2T(33.63I  j 

COMMON  /check/  ICNT.lKNfilTER 

CC  ••  change  INPUT  TO  NAMELIST 

hamelist^InUaTa/  ipBinT.ifirst.itape.tstop.capu.pscht.psitst, 

C EPS, Ra, ICNT, DELCU, REDLIN.KONTF, lasts, FACT ,F ACT l.FACT2_j 

C FtfPX,FtTP'Y,RUNNO 
CC  *• 

DIMENSION  STOREmiT^im 

C DIMENSION  US  132, 62,2)  .VS  (32. 62 ,2) J 

" ' DiHENsroN  usd.i.i) .vsa.i.i)  j 

EQUIVALENCElOT.STOREll) ) 

CC  **  ^ 

REAL  MVJO 

CC  •• 

DIMENSION  FV0(63) ,FUB(63) .FVBX(63) ,FUBX(63) ,Z(63) 

DIMENSION  An62(2l ,NAR6(2) 

C 

C PSI  IS  THE  STREAMFUNCrrON 

C U=0(PSn/0Y  » X VELOCITIES 

C V = 0(PSir/OX-»^  VECtJCITIES 

C VORTaVORTICITY 

ir ORIGTW  10. .0.)  IS  A1  l=2-,ja2 

C 

CC  *•  FUNCTIONAL^STATEHENT-rOTrEXTTTST'OCTTrOM  ■Dr  GREEN'»S  FN, 

CC  •* 


EXTRAP  (XI  iXZTn?TrTT3T~g~yiTrX2^yn7-TTe=TrT^TT3^I 
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lYM2xIY-2 

*«VJ0»32. 

CAPUSTxO.O 

CC  •*  READ  IMPUT 

READ(5tSQ03) 

5003  FORMAT UX" 

WR I TE  < 6 » SOOjn 

REA0(5«tN0ATA| 

CC  *•  WRITE  NAMELIST 

WRITE (6< INOATAI 

C PSITSTxRELAXATlON  CONVERGENCE  TEST  FOR  PSl  ITERATION 

CC  **  PSITST  SET  FOR  INITIAL  ADJUSTMENT  ONLY.  PROGRAM  THEN 

PI*3,1A15926 

RKA>3.S3r7b6 

RJOKA  — .A027S94 

PIH*Pi/'2. 

IPRTxQ 

Hal.O 

C RA  > RADIUS  OF  CIRCLE  OF  ZERO  PSI-LINE 

RKxRKA/RA 
PHAxRA 

PSCxRKA»RkA/ (flMA»RHA) 

PIINxl./PI 

'*PIITM2»-2.»PIIN 

WRITE (IPRTERt 1002)  IF IRST . I T APE .RA , TSTOP .CAPO.PSCHT tPS I T- 

lEPS 

1002  F0RMAT(1H1."THE  FOLLOWING  IDENTIFIES  PARAMETERS  FOR  THIS  RUN  - 
l//»"  IFIRST  ■ ".IZ."  WHERE«0  means  1ST  RUN".//."  ITAPE  « ".12./. 

1 RA  ■".E16.8.//."  TSTOP  » ".EI6.8.//."  CAPO  ■ ".£16.8.//."  PSCHT 

2 ".E16.8.//."  PSITST  » ".E16,8.//."  EPS  ■ ".£16.8) 

WRITE (6.6005)  ICNT 

6005  PO«H*rP' WAX.  NO.  OF  TAPE  HEC0HD5 UBITTEN  THIS  RUN — ".15/) 

A2>RA*RA 

' u2»2.*tAPU 

01»U2/RJ0KA 

C — ITape  > 0 - pAOORam  STAATS  AT'T^I)  AND  0025  NOT  WRITE  FINAL  ARRAYS — 

C ITAPE  « 1 - PROGRAM  STARTS  AT  TxO  AND  DOES  WRITE  FINAL  ARRAYS 

C FoB  TiSe^AS  INPUT  " IN  "SUBSLaUENT  RUNS 

C ITAPE  « 2 - PROGRAM  USES  PREVIOUSLY  WHITTEN  FINAL  ARRAYS  AS  STARTING  f 

C ANCwAITE?  PINAL  ARRAYS  FOR  USE  AS  INPUT  IN  SUBSEITOETIT  RUN'S 

C ITAPE  ■ 3 - PROGRAM  USES  PREVIOUSLY  WRITTEN  FINAL  ARRAYS  AS  STARTING  POINT 
C Atvo' OO^S  NOf  WRITE  PINAL  ARRAyS  FOR  uSE  AS  INPUT  Ih  PuTuRC  RuNS 

ITYPE»ITAPE*I 

C 1PRINT~VALUE  OF  PRINT  THEdDENCY 

C IFIRST  *0  FOR  FIRST  RUN 

C ■ '.NE.O  FOR  SUBSEQUENT  RUNS 

C ITAPE  * TAPE  OPTION  ( SAME  AS  RADIATIVE  BUOYANT  CONVECTION) 

C TSTOP*  TIME  TO  STOP  RUN 

C CAPU*INPUT 

CC  ••  iKNT  COUNTS  THE  NUMBER  oF  CALLS  TO  PHTOUT 

CC  •*  ICNT  « MAX.  NO.  OF  TAPE  RECORDS  TO  BE  WRITTEN  BEFORE  STOPPING 
CC  ••  'kONT  » <:OUNTS  NO.  OF  TIME  STEPS 

CC  »*  DELCU  A CALCULATED  CORRECTION  TO  CAPU 

CC  ••  OELU  - AN  ARRAY  FOR  CORRECTING  CAPU 
CC  •*  I6RNKT  GREEN«S  FN  COUNTER 

CC  ••  I6RNM~GReEN"S  PN  INOICATSR  WHICH  C6^(TAIN5  T>iE  VALUE  TJF^reRRlTI 

CC  ••  THAT  WAS  LAST  USED  WHEN  THE  FN  WAS  CALCULATED 

CC  •*  US  AN(j'vS  ARE  ARRAYS  FOR  STORING  BOUNDARY  U"S'  AnB“V""S  FOR  EX TCa'  JL**  f ION 

IF (IFIRST. NE.O)GO  TO  10 

OX*l ./60. 

OY*l ./60. 

GO  TO  n ■ ■ 

10  CONTINUE 

CC  •*  NOTE  CHANGE  TO  IFILEl 
IFILE1*10 


") 


ITST 


1 


BEST  AVAIUBLE  COPY 


I 
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[ ? IYP1*IY»1 

[ IYM2»lY-2 

I mvJ0=32. 

[ CAPUSTsO.O 

CC  ••  READ  INPUT 

REAO(St5003) 

1 5003  FORMAT  UX”  ■<! 

i WRITE(6.500JLI 

REA0(5.IN0ATA) 

CC  ••  WRITE  NAMELIST 

WRITE (6tIN0ATA) 

C_  PSITST»RELAXATI0N  CONVERSENCE  test  for  PSl  ITERATION 

I CC  •*  PSITST  SET  FOR  INITIAL  ADJUSTMENT  ONLY,  PROGRAM  THEN  RESETS  PSITST 

PI»3.1Ai5926 

RKAa3.83r7b'6 

S RJOKA*-, *027594 

! PIH»Pi/2V 

IPRT»0 

H*  I , 0 

C RA  « RADIUS  OF  CIRCLE  OF  ZERO  PSI-LINE 

RKzRKA/RA 
PHAsRA 

PSC»RKAiRkA/<rtHA»RHAI 
PI IN*1./PI 

[ PIITM23-2.*PIIN  ^ 

S WRITE(IPRTERfl002)  IF IRST 1 1 TAPE tRA , TSTOP «CAPO, PSCHT.PS I TST , 

lEPS  ‘ ' ~ ■ ^ 

1002  FORMAT (IHI, "THE  FOLLOWING  IDENTIFIES  PARAMETERS  FOR  THIS  RUN  - ■«, 

[ l//f“  IFIRST  » '•,12.”  WHEHEaO  MEANS  1ST  RUN",//,"  ITAPE  » ",I2,//," 

r I RA  3",EI6.0,//,"  TSTOP  » ",E16.B,//,"  CAPU  « ",E16.8,//,"  PSCHT  » 

I 2 ",E16.8,//,"  PSITST  » ”,£16.8,//,"  EPS  « ",E16.8) 

WRITE (6,60051  ICNT 

6005  FoRMArn'  MAX'.  NO.  OF  TAPE  RECORDS  WRITTEN  THIS  RUN — R;I5/) 

A2sRA»RA 

*CApu  ~ “ 

01SU2/RJ0KA 

C ITAPE  a 0 - PROGRAM  STARTS  AT  T»tf"JWD  ODES  WT  WRTTE  FINAL  APRAVS 

C ITAPE  " 1 - PROGRAM  STARTS  AT  TaO  AND  DOES  WRITE  FINAL  ARRAYS 

C FoiT  uSe  AS  INPUT  IN  SUUSEOUENT  RUNS 

C ITAPE  a 2 - PROGRAM  USES  PREVIOUSLY  WRITTEN  FINAL  ARRAYS  AS  STARTING  POINT 

C an6 "writes  FTMAL  ARRAYS  FOR  05r~AS  INPUT~rTr  5CrH5Eni3ENr“R0NS 

C ITAPE  a 3 - PROGRAM  USES  PREVIOUSLY  WRITTEN  FINAL  ARRAYS  AS  STARTING  POINT 

C AND  DOES  NOt  WRITE  FINAL  ARRAYS  FOR  u5E  AS  INPOT  IW FuTo1»E' ROWS” 

ITYPEalTAPE*! 

C 1PRINT~VALUE  OF  PRINI  FHEQUENCT 

C IFIRST  ao  FOR  FIRST  RUN 

C ".NE.O  FOR  SUBSEQUENT  RUNS 

C ITAPE  a tape  OPTION  ( SAME  AS  RADIATIVE  BUOYANT  CONVECTION) 

C TSTOPa  TIME  TO  STOP  RUN 

C CAPUaINPUT 

CC  ••  IKNT  COUNTS  THE  NUMBER  OF  CALLS  TO  PRTOUT 

CC  •*  ICNT  a max.  no.  OF  TAPE  RECORDS  TO  BE  WRITTEN  BEFORE  STOPPING 
CC  *•  kONT  a COUNTS  NO.  OF  TIME  STEPS 

CC »*  OELCU  A CALCULATED  CORRECTION  TO  CAPU 

CC  •*  OELU  - AN  ARRAY  FOR  CORRECTING  CAPU 
CC  •*  IGRNKT  GREEN”S  FN  COUNTER 

CC  lGRNH~MEgWi'g'  F><  I^tOlCATOft  WHICH  CONTAINS  TRf  VALUE  OF~rCffRirT 

CC  ••  THAT  WAS  LAST  USED  WHEN  THE  FN  WAS  CALCULATED 

CC  •*  US  AND  VS  ARE  ARRAYS  FOR  STORING  BOUNDARY  U”S  ANC"  V""S  F^  EXTRAPOLATION 

IF(1F1RST.NE.O)GO  TO  10 

DXal./60. 

OYal  ./60. 

GO  TO  11  ■ ■' 

10  CONTINUE 

CC  ••  NOTE  C)4aNGE  to  IFILEI 
IFlLElalO 


BEST  AVAIWBLE  COPY 
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••  CHANGE  TO  GET  THE  LAST  "STORE"  OF  THOSE  "STORES"  PREVIOUSLY  WRITTEN 

**1-20 

**  2 - *0 

**  3-60 

**  A - 80 

••  5 -naa 

**  6 - 120 

**  7 - lAO  ■ ■“  ~ 

**  8 - 160 
**  9 - 180 

**  10  - 200 

DO  15  r«l  .LISTS 

REAOIIFILEII  STORE 
15  CONTINUE  ' 

•• 

11  0X1  = 1. /OX 

oyi=i./oy 

0Y2A2;*DY 

0X2=2. *Ox 

0XIH«. 5*0X1  

DYIH».5*0YI 

0X21=0X1*0X1 

0Y2I*0YI*0YI 

- 0R2»0X*0X*0T*0T 

0R2I=1./0R2 

T0=0X2I*0T2T— 

TD2»2.*T0 

FMU» I0X2I/T0) *0057? I "OX  J * (0Y2I /T  01 'COS  IP  I •OY) 

A » FHU/n.*SORT(l.-FHU*FHU)  1 

■ PL*lT-rl.«A*A» 

PC=(1.0*A*A)*.5*OR21 

P1=0Y*0Y»PC 

P2=DX*0X*PC 

P3*0T*0Y*P2 ^ 

•• 

TrrTFTHST.Nt.OJ  00  lU  99 

00  20  I=I,IXP1 

DO  20  JAITIYPI — 

U(IfJ)=0.0 
VIIiJ)*0.'0 
PSI (I«J) =0.0 

voRT(T*.nTi»o;o 

VORTd.J, 21=0.0 

20  CONTINUE  

Cl  = (2^0*C*PUJ/{RK»RJ010n 

TIME=0. 

0T*0.  

SET  UP  UPPER -OlIAOPANT  OF-rmTI*l.~Y»CUrS“T)F— CtOSElT  CIRCLE  TWO 

.USE  SYHMETRY  for  BOTTOM  QUADRANT 

II»2 

I2»IIX-2|*. 5*2. 0009 

Jl  = (lY-2)*. 5*270009 

J2«(IY-2)*. 75*2. 0009 

J3«(IY-2)*. 25*270009- — 

00  30  1=11.12 

00  30  J»jriJ2 

X»FL0AT (1-2) *0X 

Y=FL0AT(J-2)*0Y-i5 --  

R=SORT IX*X*Y*Y) 

PXR=R*R  — 

SINTM=X/R 


BEST  AVAILABLE  COPY 


V 
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toSlHa-r/H 
IF (R-RAI 23,23,25 
23  ARG=Rk*R 

CALL  BESJ(AR6,1,BJ,. 00001, lERI 

I F ( I CR . Nrr5TWRlteiT;5099JTtR,I,JaR6,SJ 


4R99  FORMATdH  ,''1ER  = ''.la,"  I*  . 1 2 , J»  ",I2,"  ARG»  »,E16.8, 


i '•  BJ  » ••,E16.B) 

PSI (I,J)»-CI»BJ*SINTH 

VR=-§I» (BJ^AHG)«COSYh  ^ 

CALL  BESJ(AR6,0,BJ0,. 00001, lER) 

IF (ftR.NE. 0) WRITE (4,4989 1 IER,I,Jt ARG,S JO 

4989 

FORMATdH  ,”IER  a '•,12,”  la  ”,I2,”  J«  ",I2,"  AR6  ■ ” 

,E16.B, 

1 » BJd  a ••,E16.e) 

VTHETa6l*SINTH»(BJ0-BJ/ARG) 

U ( i , Ji  •♦VRaS  INTH,VTHE:T«C<JSTH 

V d , J) a-VR*COSTH,WTHET*SlNTH 

VbRt ( I , J, rrSPS^aRSl (!,J1 

V0RTd,J,2)aV0RTd,J,l) 

60  T()  26 

25 

PSI d,J|a-CAPU*(R-A2/R)*SlNTM 

VRa-CAPU* ( 1 .-A2/RXR)»C0STH 

VTHETaCAPU* ( 1 . ♦A2/RXR) *SINTM 

Ud,Jia»VR*SINTH»VTRET*COSTH 

W (1 , J) a-VR*COSTH,VTHET*S  INTH 

26 

KaJl-(J-Jl) 

PSId,K)aPSI  d,J) 

Ud,K)a-u<I,J| 

Vd,K)aWd,J) 

VORTd,K,l)aVORTd,J,l) 

V0RTd,K.2)aV0RTd,K,ll 

30 

CONTINUE 

C 

INSURE  THAT  GRID  POINTS  ON  CLOSED  CIRCULAR  BOUNDARY 

ARE  ZERO 

WRitE(6,4979)PSI  (n,Jli  ,R4l  d2,Jll  ,PSI  dl,J2>  ,PSI  (iT 

,J3) 

4979 

FORMATdH  ,4E16.8) 

PSI  (Il,JI>aO.O 

PSI  d2,Jl)a0.0 

PSIdl,J2)aO.O 

PSI(Il,J3)a0.0 

C 

SET  UP  REST  OF  INITIAL  PSI  FIELD 

I3aI2*I 

DO  40  IaI3,IXPl 

00  40  Jajl,IYPl 

X a FL0ATd-2l*0X 

YaFL0ATtJ-2)*0Y-.5 

R a SORt(X*X*Y*Y> 

SINTH  aX/R 

COSTMa-Y/R 

RXRaR»R 

PSI d,J>a-CAPU*(P-A2/R)*SINTH 

VRa-CAPU* ( 1 .-A2/RXR) aCOSTH 

VTHETaCARu*d.*A2/'RXR>  *5 INTH 
Ud,J»a*VR*SINTH*VTHET*COSTH 

V ( 1 ,J>a-vR*COSTH»VTHET»SINTH 

KaJl-(J-Jl) 

PSI (t,K>aPSI (I,J> 

Vd,K)aVd,J) 

40 

iMNKla-oriVJT 

CONTINUE 

J4aj2»l 

00  SO  Iall,I2 

00  40  Jaj4,tYPl 

XaFL0AT(t-2)*0X 

YaFL0ATfJ-2)»DV-."5" 

WaSORT(X*X*Y*YI 

PXRaPWR 

SINTHaX/R 

BEST  AVAIUBLE  COPY 


r 
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t 


f 


I 


CosTM=-y'/(^ 

PS  I ( 1 1 J»  =-CAP0*(P-A2/P| 'SINTH 
VR=-CAPU* ( l.-AZ/PXR) “COSTH 
VTHET*CAPU» ( 1 .♦A2/RXR) ‘SINTH 
U(t.JI*»VR*SINTH«VTHET»COSTH  " 

V ( I t J) »-VR«COSTH*VTHET»SINTH 

K=Jl-(J--jll  — 

U<lfKI*~u<ItJI 

V(I.K)aV(ItJ) 

PSI (l.K>»PSI (ItJ) 

50  CONTINUE  ' - 

CC  ••  constants  defined  on  the  boundary  BG  (BEFORE  GREEN) 

c psiTP*P5n37Tin 

C PSIRTaPSI (30,32) 

CC  *•  ■ " ■ 

00  GO  J=1,IYP1 
GO  PSI (l,J)a-PSti3,J) 

DO  TO  J=1,IYP1 

VOPT ( 1 ,J,l>a-V0Rt(3,J,l) 

V0RT(l,J,2)aV0RT(l,J,l) 

TO  CONTINUE  " 

C CALCULATE  U,V  FROM  PSI  FIELD 

C PSI  IS  THE  STREAHFONCrrON 

C UMAXtVMAX  ARE  USED  TO  ADJUST  DT 

UMAjfiTri 
VMAX*2.5 
UMAXa2A. 

VMAXaSO. 

UMAX«1200. 

VMAXa2S00. 

CC  *•  RESET  P5ITST(PSI  TESII  FOR  U5E~IM  FIRST  STEP  0NL71 

IF(  KONT.EO.O)  PSITSTaPSITST*.! 

CC  ••  

CC  •• 

CC  ••  ” 

CC  •*  HERE  IF  STARTING  FROM  TAPE 

9<r-C0NTIITOE 

MRITE(G,INOATA) 

CC  •*  GET  INITIAC-OETAICEtr-PRINTOOT 

C CALL  PRTOUT 

CC  *•  

CC  ** 

-100  CONTINUE 

CC 

C 

CC  •* 

230  0TIX=UMAX«OYI»VMAX«OX1*TO2 — 

0TCRIT«1./0TIX 

OTaiB'OTCR-IT 

CC  *•  modify  OT  FOR  first  STEP  ONLY 

IF(KONT<EO«0)-  l>T»OT*TOi 

C IF(KONT.EO.O)  OT»OT».l 

CC  **  — 

240  TIME*TIME»0T 

•RITt(6,232») 

2325  FORMAT)///) 

WRITE (6,2300)UMAX^VMAX 

2300  FORMATdH  ."UMAX  * ".EIG.B,"  YMAX  a ".EIG.B) 

WRITE(G,230UDTrTIME- 

2301  FORMATdH  ,"0T  a ", EIG.B,"  TIME  a ", EIG.B) 

CC  ••  

CC  ••  CAROS  SPECIFING  BOUNOARY  VORTICITY  values  go  HERE 

CC  ••  — — — - — 

CC  ••  CALCULTAE  NEW  VORTICHIES  EVERYWHERE.  EXCEPT  AT  EXTERIOR  POINTS 

00  190  Ja2,tY  

00  190  Ia3,IX 


BESl  AVAlWBlt  cm 
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IF  (V1.6T.O.OIGO  TO  180 


V2sVl*V0RT(I»l.Jtn 
eO  TO  181 


180  V2-Vl*VORT(IiJtl) 

181  V3«U(I.J)»U(I-1»J1 


T5r 


182 

183 


IF  IV3.6T. 0.0)60 

V*»V3«VORT lltj.l) 

60  TO  18‘3 

V»»V3»V0RT(l-l,J.l) 


VS«VI1*J*1)*VII.J) 
IF (V5.0T. 0.0)60  TO 


l&L. 


V6«V5*V0RT(I»J*1*1I 
GO  TO  185_ 


18*  V6aV5»V0RTII«JtI) 
185  V7«VI1>J)*VII>J-1) 


IF(V7.6T. 0.0)00  TO  186 
V8-V7»V0RT(ltJtl) 


60  TO  187 

186  V8»V7»V0RTI1«J-1.1) 


187  VClaDXlH*(V2*V6) 
VC2a0YIH*«V6-V8) 


VC1T(I.J)«VCI 

VC2T<I.J)»VC2 


VC«2.*V0RT(ItJ.l) 

VC6«0X2I*(V0RT (I*l.J.l)-VC*VORTII-l.J.l) ) 
VCS«dY?I*(V6Rt(ItJ*l.l)-V(i*V0HT(l,J-l,l) t 


CC  •• 

cc  •• 

CC  •• 
CC  •• 


SIONU  - INOICXTES  WMXT  VALUE  TO  USE  IN  CORRECTING  FOR  FXLSE  VISCOSITY 


• SI6NV  - INDICATES  WHAT  VALUE  TO  USE 
SI0NU>SI6N(.$t-U(If J) ) 
SI6NV-Sf6N(.5fV(NJ)  ) 

IFIU(l.J).EQ.O)  SIGNUaQ. 

IF  fw‘ ( iVJ)  .EQ.O)  $I6NV«0. 


In  CORREdTINO  FOR  ^ALSE  VISCOSITY 


CC** 


TvCIl ^iSNU*D)(«(u(I.J)*VC4*TPC01(0n*l«Jt.U(I-l.J).0Xi  * 

: TPC01<V0RT(I*l.J.l) .VORT(I-l,J.l) .OX) ) 

FVCYisi5NV*CT* (V(I.J)*VC5»TPC01(V(I.J*IT7VTI.J-1).0Y)» 

TPCOl  (VORTtI.J*l»l)  .VORTU.J-l.l)  .DY) ) 


CC 

CC 


C 

• • 
• • 


FITPX  ANO  FITPY  ARE  FITTING  PARAMETERS  RELATED  TO  THE 

PALSn)lSrt5IT7"C0RRECrnJN  TEBHS  IR  X KUO  IN  V — ANO'TWE 

NORMALLY  SET  EQUAL  TO  1 


CC  •• 
CC  *• 
CC  •• 

CC 


•• 


t $ s s 


FCFX(I.J)«FVCX 

FCFY(1.J)«FVCY 


IFIABSIFVCX) .6T.FITPX*ABS(VC1) ) FCFX ( 1 . 0) *FCFX t I . J) *1 .E20 
IF  « ABS (FVCY) »6T.FlTPY*ABS<yC2) ) ^ < I ♦ "^CFY ( 1 1 J) * 1 .E20 

IFlABSiFVCY) .GT.FITPY*ABS(VC2) ) FVCYaSIGN( VC2.FVCY) 

VORt ( I .^.2) -VORT ( 1 , J, I ) ♦0T»rr7C**VC5-VCr-VC7T 

NOTE  V»S.0U1  .\RM  lE.e.OAtS)?  8 ^ 

V0RT(I«J.2)aV0RT (I .J.l  ) ‘OT* ( < 0 , 0»0 . 0-VC 1 -VC2 ) 


C 

..8  ** 


♦FVCY) 


C ‘FVCX 

190  CONTINUE' 

00  107  I«3.IX  

V0RT(1.1YP1»2)»3.*(V0PT)1»1YP1-1 t2>-V0RT(l.lYPl-2.2) ) 


C ♦VORT(I.IYPl-3.2) 

CC  ••  CALCULATE  VORTICITY  ALONG  BOTTOM  EXTERIOR  LINE 

VORT  I i . i . 2 ) i3  . * 1 VORT  i r. 2 .'2 ) - VOHt  ( 1 . 3.2 ) ) • VOfi T ( 1 i * » 21 
107  CONTINUE 

CC  ••  calculate  VORTiCITY  ALONG  H ITjHTTxTER lOR  LINC 
no  108  J«1.IYP1 


BEST  AVAIUBIE  COPY 
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VORT (IXPl ,J,2) *3.»(V0RT (IXPl-ltJ,2)-V0RT  nXPl-atJ,2) ) 
C VOPT(IXPl-3fJt2) 

108  CONTINUE  ■■ 

CC  •• 

DO  1<J1  J»liIYPl  

VORT(ltJt2)>-VOPT(3.J>2) 

191  CONTINUE  ■ 

CC  •*  BY-PXSS  eREEN"S  FN,  EXTRAPOLATION 

CC  •*  t(HEN  STARTING  FROM  TIHrTERO  - _ . 


IFIKONT.LE.Z0)  GO  TO  195 


CC  •*  TO  AVOID  SKIPPING 

IFIKONT.GE.O)  GO  TO  195 

CC  ••  

CC  •*  lORNKT  IS  COUNT  CONTROL  FOR  CALCULATION 

CC  *•  OF-OREEW"S  FUNCTTOW 

CC  *• 

CC  •*  always  PERFOR)r‘CACCUC)lTTOI9'WHEN  STARTING  OFF 

CC  •* 

IGRNKTalGRNKTXl — — — 

CC  ** 

- — IFfHODdORNKI.J)  .EU.0.ANU.UT.Lr;r.5L-^rGP~T0~575 

195  CONTINUE 

WRITE IA«601OI 

6010  FORMAT!"  GREENS  FUNCTION  USED  DURING  THIS  STEP  ") 

C IMPOSE  OREENS-FONCTrONS-Ttr-OrrAIN-UTV“WCOONOARIEy 

C 

C - SEABCfT  FDR  NON-ZERO  VORTTCITIES 

C 

NVU»0  

NYL«0 

NXR»0  ^ 

NXL«3 

- DO  9G-a*JTTTMT ^ 

IF(ABS(V0RT(3.J.l) I.LE.EPS)60  TO  90 

NYL»J  

60  TO  91 

90  CONTINUE  

91  DO  92  J»1,IYM1 

K5IT=-J 

IF  (ABS<V0RT)3tK«l) ) .LE.EPSI60  TO  92 

NYU»K  

GO  TO  93 

92  CONTINUE  ■"  — “ 

93  DO  96  J»NYL»NYU 

DO  94  lijmmi 

IF(ABS(V0RT(I.J.1>).6T.EPS)00  TO  94 

NX*I*1 — 

60  TO  95 

94  CONTINUE  ~ 

95  IF (NX.GT.NXR)NXRsNX 

96  CONTINUr  

WRITE ( IPRTER. 1050) NXL.NXR.NYL. NYU 

1050  FORMAT(lHOf"  NXir»“"Tr2T"‘NX(r-w-"Tl27"^TL'»  ".12."  NYU  « ".12J 

C 

C LIMITS  OF  NON-ZERO  VORTICTTr-AREA-MAVE  BEEN  ESTABLISHED  " 

C 

C CALCULATE'U.V  ON  Y*(J“)INU  T«1  BOUNUARV  LINES 

C 

ABG2ri)a0.0  - 

APG2(2)aIY-2 

NARG(l)a2  - • 

NAR6(2)sIY 


BtSl  AVAllABlt  COPY 
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00  5»0  L»1»2 
YB»ARe2(L» 


C 

cc 

cc_ 

cc 

cc 

cc 

cc 


t.L*N*p5(Ll 

D0^500^n«3t 


JSkfP  . \ 6k  i be«kiOIN6  wHETWEB  6(»EE^<'•5  TN.  15  INTE RP0LATE6 
AT  EVERT  OTHER  POINT  OR  CALCULATED  AT  EVERT  POINT 


JSKIPaZ 


00  500  lla2«lX«JSKIP 


K = l 

X8aFL0AT(n-2» 

xe2>XB*XB 

FVBX(n*0.0 


FUBX<I)«OTO 
DO  *20  J»NTL»NTU 


KxiTfT 
YP«FL0AT(J-2) 

FVB(n=0.0 

FUB(1)30.0 

FV8tNXR)»0.0- 

FUB(NXR)aO.O 


CY»YB-YP 
CY2aCY*CT 
Ob  A 1 0 T¥3 . NXR 
XP»Ft.0AT<t-2l 

OENbMV(  (X8^-XP>»(XB-XP)  ♦CYE  )»<  IXe*XP»  * (Xfl*XP)  ♦CyJ  T 
FVBC1-1)»V0RT(1,J.1)»0XI»XPA((XB2  -XP*XP-CT2  )/OENOH) 


AlO 


FOB  ( 1-1 ) r ( VORT  < 1 1 J.  1 » *6xi*xa»xJ‘*CY)  /OENOM 

CONTINUE  

CALL  OSPToXfFVBtZ.Tixfn 
FVBX(K)»Z(NXR) 


CALL  QSF(OXtFUB«ZtNXRI 
FUBX(K)»Z(NXH1 


420  CONTINUE 

K*K*1  

FVBXfKjioT? 

FUHX  fKJf  0_.  0 

CALL  QSFjOTtFVBXfZ.K) 
VllltLL)«Z(K)»PllN»CAPU 


500 


CALL  OSF (OTfFUBXtZtK) 
U ( 1 1 tLL)aZ(K)»PllTH2 
CONTINUE 


CC 

CC 


FIRST  SET  OF 
FIRSt  SET 


U ANO  V MOOIFIERSGO  HERE 


00  SOS  Il33tIXMlt2 

U ( n .LLJ  a.( Um-I  «LL)*U(  11*1  «LL>>  4.5 
505  V(n«LL>>rvin-ltLLi*V(Il*ltLLI  >4.5 

CC  4* 

510  conTInoI  " 

C 


calculate  u*v  ON  x».5  LinE 


APG2(Uai) -0"~ 
APG2(2)>IX-2 

nargh»*2  — 

NARG(2)»IX 


LL»Z 
II»NAPG(LL» 

X8*AP62Ill1 

DO  550  L«3« lYMl.JSKIP 
K*1  - - 

FVBXtl)»0.0 


%©  m 


j 


f 


i 

i 

i 

( 


i 

I 

I 

i 

! 


t 


t 
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ruBx 1 1 » *0.0 

Ye*FLOAT (L-2) 

00  540  J»NYLtNYU  : " ‘ 

KaK*l 

YPsFLOAT(J-H)  - “ 

CY»Ye-YP 

CY2«CT«CY 

00  530  I*3«NXR 

XP«FL0AT(l-2)  

DCN0M«( (XB-XP)*tXB-XP) ♦CY2  > • < (XB*XP ) • (XB*XP) *CY2  ) 

FV8 ( i-i I svoRTrn jrn *(m*xp»TrxHZ — =xp»xp=cyz  ^ r/ozNOMj 

FU8 ( I-l ) « ( VORT ( I t Jt 1 ) *0X l*XB*XP»CY ) /OENOH 

53(r  COMTINOE 

CALL  OSF(DXfFVB<ZfNXR) 

FVBX (K » *Z  INXR  J 

CALL  OSF(OX.FUBfZ*NXR) 

FUBXCKJsZINXHJ 
540  CONTINUE 

■ fTifcri 

FVBX (K)«0.0 
FUBXlKIsr.O 
CALL  OSF(OY.FVBXtZ.K) 

V(IliL)i2(KI»PITWiC7IPO 
CALL  0SF(DYtFU8X.ZtK) 

UTTriC>*Z<IO»PnTH2 

550  CONTINUE 

cc  •* 

CC  ••  GREENS  FUNCTION  TEST 

lYN2«lYHr-l 
CC  *•  SECOND  SET  OF 

CC  *»^  (riNO  V M60inc«s  06  hErE 

CC  **  SECOND  SET 

00  555  C**tITH272 
U(IXtL)>(UCIXtL-l>  *U(IXtL*l) )*.5 

V ( I X .L I « { V ( rXTC=lTFVnXiL  * rTT*75 

555  CONTINUE 


CC  

60  TO  (5S8«588) tIONE 


CC 

c 

CC 

575  CONTINUE 

• # 

CC 

«• 

GO  TO  (590#590) tIONE 

CC 

CC 

• • 

CC 

CC 

• * CHANGE  V (AFTER  I TO  V tBEFOREI  

• • 

NDIF*IT-2 

00  580  I«2tIX 

00  580  J«2t62fNOIF 
V(ItJ|aV(ItU>-DELCU 

^80  CONTINUE  ' --  --- 

00  se5  J*3tIYMl 

v(iXtJi*vrtXfUr=OET:CT) 

585  CONTINUE 

CC  •*  END  CHANGE  “ 

588  CONTINUE 

590  CONTINUE  • - - 

CC  ** 

C calculate  PSl-ON  BOUNDIWIES-  AT  TO^T80rTOW-ANO~AT~ffTqwT 
C 

PSI  (3.2)«-0X*V  (2*2> 

PSI (3.IY)s-DX*V(2tIY) 

00  700  I»3*IXM1  - 

PSI(I*lf IY)=PSI(I-l.lY)-0X2*Vtt*tY) 


aCi 


BEST  AVAILABLE  COPY 

J 
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PSI(I»1.2)=PSI(I-I.2)-DX2*V(I.2) 

700  CONTINUE 

CC 

PS  I ( I X . 1 ymi)=2.*psi  ( ix.iY)-psi  uxiTitiv)  *ox»vnTrrTr=DY»urix.rYJ 

•*  CHANGE  TO  REFLECT  VORTICITY  ON  THE  BOUNDARY 

C -.5*OX*OX*VOPTMx.rY.2J 

00  710  J»3.IYM1 

k»1ypi-j 

PSHIX.K)»PSI  <IX.K*2>-0Y2»U<IX.K*1) 

710  CONTINUE 

250  ITER«0 

260  DPSIMXxO. 

PSIMAX«0. 

CC 

ITER-ITER*! 

**  TEST  ITER  TO  PREVENT  RUN  AWAY  ITERATIONS 

IF(ITER.GE.35)  CALL  PRTOUT 

00  270  J>3.IYM1 

00  270  Is3.tXMl 

PSI (I.J)  »RL*PSI (I.J) ♦P1*(PSI (I*1.J»*PSI (I-l.J) )*P2*(PSI  (I.J*1) ♦ 

IPSI (I.J-1) )»V0RTCI.J.2)*P3 

270  CONTINUE 

DO  280  Ix3.IXM1 

DO  2B0  J>3tIYMl 

PS«PS1 (I,J» 

PSI (I tJI  *RL»PSI (I.J) ‘Pl^IPSI (I*1.J)*PSI (I-l.J) )*P2*(PSI (I»J*1) * 

IPSI (I.J-1) )*V0RT(I.J.2)*P3 

OPSI>ABS(PS-PSI (I.J) ) 

CC 

IF(OPS1«LE.OPSIMX)60  TO  288 

CC 

OPSiMXaOPSl 

••  FIND  coordinates  OF  OPSIHX 

CC 

I Turn 
nH2*j 

A*  END  FINO 

288  CONTINUE 

( ABS  (PS  I ( 1 9 J H ) 60  10  <!tlO 

PSlHAXsABSCPStdfJ)) 

CC 

280  CONTINUE 

•*  CHANGE  J BOUND  TO  GET  CORNER  POINT 

00  290  j«irrv 

PSI (I.J)x-PSI (3.J) 

CC 

290  CONTINDE 

•*  PSITST  TEST 

CC 

WRlTE(69600n  IIHlf  IIHZ 

6001  FORHArrAOX?2TTir) 

WRITE(6.6002)  OPSIHXfPSIHAX 

6002  FdRM/iri^F^i) 

CC  *• 

c 

irTOPSIMx.GT.PSITSTAPSIMAX)GO  10  260 

CC 

♦A  RESET  Psi  test  VALUt  i-oh  all  SUHStOUtN)  FINt  STEPS 

IF(KONT.EO.O)  PSITST»PSITST*10. 

CC 

'•c 

• # 

**  WRITE  INOATA  REFLECTING  CHANGE  OF  PSI  TEST 

CC 

If  (KONT  *£0  • 0 ) 1 1 E ( 6*  INUA  1 A ) 

• « 

DO  281  Ixl.IXPl 

DO  281  Jxl.IYPl 

V0RT(Iij;i7SV0RT(T.J12T~  

281  CONTINUE 

c 

CC  •*  COUNT  NO.  OF  TIME  STEPS 


K0NT»K0NT*1  ' * ~ ■ 

MRITEI6.6022)  KONT 

6022  FOHM*t(w  — PJUWBEir"(5F' TWrS  'TTME  S7CIF  ."  -".IST 

C OETEPMINF  change  in  CAPU  * CHANGE  IN  V OF  STAGNATION  POINT 

BtSl  AVAllABlt  COPY 


1 

I 


I 


i 

i 
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1,1,  ••  l<^T  MAX  PSD  O1VI0E.0  BY  DT 

CC  •*  _ _ 

C 

MVIaO  _ 

MVJ»0 

PSIMX«0. 

00  600  IV37T)( 

00  600  Ja3iIY 

IF(PS1  (I.J)  .LE.PSlMxrSO  TO  600 
MVI»I 

MVJaJ  

PSIHXaPSKDJ) 

600  CONTINUE  ■ 

CC  **  FIND  max  PSl  IN  THE  NEIGHBORHOOD  OF  PSKMVItMVJI 

CC  **  — 

CC  ••  USE  THE  BEST  NINE  POINTS  ANO  THE  SECOND  BEST  NINE  POINTS 

CC  ••  ■ 

CC  •* 

~ MD*SIGNT1  . « (PSITWVI  ♦ I f HVJ)~PS1  lHVI-1  ,HVJ»  > » 

KD>0 

610  CONTINUr 

A0I»PSI  (MVI.MVJ*n-PSI  (MVI.MVJ) 

AOMIaPSr(HVr*MVU-D-PSl  IMVI.HVJ) 

A10>PSI (HVI*lfMVJ)-PSI (HVItMVJ) 

-AHrOaPSI  (HVl-l.HVJI-PST  IHVnHVJl 

CC 

B5«.25»«PSlCMVTFITMVUm 

1 -PSI IMVl-l.MVJ*!) 

1 - -Psr(>iviTrni7-j=Ti 

1 ♦PSI tMVI-l,MVJ-l> ) 

Bir;5»TAn»-AH10» 

B2».5*<A01-A0M1) 

83*i5*(AlO»A>Yra) 

BA».5*«A01*A0MU 

0ELTA»Ai*B3»B*=B5*B5 

0ELX»-t2.*Bl»BA-B2*B5)/0ELTA 

0ECT^t'‘B2-B5»DeLX»».3/8A 

DELPSI  =B1»0ELX*B2*0ELY^B3*0ELX*»2*BA«0ELY«*2»B5*0ELX»0ELY 

PSIMXaPSr(MVDMTUr»OCCPSl 

WRITE! 6*6025 > PSI  (HVItMVJ)  *PSI  (MV I tMVJ*  1 ) «PSI  (MVDMVJ-1)  , 

c PSI  (Mvi»r*iFvynrPsrTMTn=T*Mvj^» 

WRITE (6*6025)  PSI (MVDl *MVJ*1) *PSI (MVI ♦! *MVJ-1 ) *PSI (MVI-1 «MV J-1 ) > 

- - C PmHVI-ltHVJ*!! 

WR1TE(6*602S)  OELXtOELYtOELPSI 

6025  FORMAT (5E20Tr) 

IF)KO.E0.1)  GO  TO  650 

KO«l  

W1»I.-ABS!0ELX) 

0ELX0iMn»DEl.T(  ' 

delyo»dely 

PSIMXO*PSIMX 

MVIsMVDMO 

iF(Ko.Ea.i)  GCTTff  mr 

650  CONTINUE 

DELxsHvnoepr  ' 

CC  *• 

CC  ••  - — 

CC  •• 

CC  •*  THESE  EXPRESSIONS  WEIGHT  THE  RESUUr  TOWARD  THE  MOST  CENTRAL  ONE 

CC  •* 

0ELXP-W1»0EL“X0T(  1 .-wn  WUEEX 

0ELX»0ELXP-MVI 

psimxspsimx*mo»oelx»(psimx-psimxo)  ■■ 

CC  ••  this  Expression  weights  the  result  toward  the  most  central  one 

OFLYPP«MVJ»OeLY*MO»OELX  • (OELY-DELTO ) 

CC  *• 


BESl  AVAIIABIE  COPY 
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cc  *• 

C MRITE (6<A02S)  AOltAOHltAlO.AMIO 

MRITE(6t6025)  AOltAOHltAlO.AHlO 
C WRITE (6f 6025)  B1 • B2 <63 t BA t B5 

wRlTE(5.6d25)  01*S?7B7T01TS5 
WRITE  15*3010)  PSIMX*DELXP.DELYPP 

3010  FORMAT  ("O"*"  MAX  PSI  a*'tE16,a,"  AT  1=  ".Pd. A."  j3  ",F0.A) 

cc  •• 
cc  ** 

cc  ••  2 POINT  PREDICTOR 

CC  ••  A PbtNf  CQrreCTQR 

1F<K0NT.GT.3)  GO  TO  620 
OELU(KONT) ■OELYPP-RtbLIN 
60  TO  6T0 
620  CONTINUE 

0ELU(I)*0ELU(2) 

0ELU(2Ja0Et.U  (T) 

OELU ( 3) =OEL YPP-REOL IN 
WR  rTF(6*6030)  (DELU(t) .1*1*3) 

6030  FORMAT!//"  YS" .5X . AF 15 .7// ) 

IF(KONT.Trr.^OI  60  T0'670 

YPRED»2*0ELU (3) -OELU (2) 
yBES=(32.*DELu(S)-12.»dElU(2) )/2t, 

cc  **  FACT  IS  A PARAMETER  TO  CONTROL  OSCILATIONS  IN  THE  CONTROL  SYS. 
CC  AS  IS  FACTl 

CC  ••  AS  IS  FACT2 

CC  »*  THE  FOLXOWING  TACTS' ARE  FOR  CAPU(O) =1000, 

C FACTs^ 

C FACTl*. 303 

C FACT2«1. 

CC  *• 

CC  ••  THE  FOLLOWING  FACTS  ARE  FOR  CAPU(0)*20, 

C FACTi7? 

C FACT1*0 

C FACf2»A 

CC  *• 

OELCUF*,5*(YOE5-YPREO)*OY/OT/FaCT2-  fact* (OELU «3) -OELU (2) I^By/CT 
C -FACT1»(DELU13)-2.»0ELU(2>*DELU(1)1»0Y/0T 

CC  *• 

CC  **  

C DELCU*OELCUC*OELCUF 

CC  *»  NOTE  ONLY  DELCUF  USED 

0ELCU*0ELCUF 

CAPU*CAPu*OELCU 

67b  CONTINUE 


WR I TEJ6 . 30 13)  CAPU 

3013  FORMAT (//"  CAPU*" .El 6. 7) 

WR1TE(6,3Q111  DELCU 

3011  FORMAT (IM  DELTA  CAPU  a"»E16.7) 

CAPUST*OELCU 

CC  *• 


1 

CC 

cc 

cc 

cc 

cc 


WRlTEdPRTER.  1008)  ITER 

008  format  {lH  .7i  ITER  * ".13  ) 

• » 

- - 


alter  PsrTo 


“(CAPosn 


DO  600  li3TI)r 

DO  580  Ja2.IY 

PSI  (1  .J)=PS1  (l.J)-tAFuST»OX*n-2) 


680  CONTINUE 

DO  602  jwa.ir 
PSt (l.J)a-PSl (3,J) 
682  CONTINUE 
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CL  “*  alter  V ON  eoUNOARY  TO  ftEFLECr  CHANGE  IN  FREE  STREAM  VELOCITY 
DO  690  I=?tIX 

V(It2)=V(lt2)*CAPUST  

V(1.IY)3V(1»IY)*CAPUST 

690  CONTINUE  ' ~ ‘ 

00  692  Js3fIYMl 

V ( 1 X t j T s V T rrrjj*cAPU5T 

692  CONTINUE 

CC  ••  END  alter  ■ ' 

CC  ** 

CC  **  U AND  V GENERATEO'BV'FCO " 

CC  •• 

- uHAx»a;o 

VHAXzO.O 
00  110  Ia2iTXHI 
DO  no  Ja3tIYMl 

u ( 1 * j I AOY I H*  iFSTt  rTjFi  T ^snitJ-m 

UTRY=A8SIUII*J) ) 

I F {OTRVTGT.UMAX>UHAX»U'rRY 

V( It J) »-OXlH»(PSI ( l*lt J)-PSI (I-lt J) ) 

VTRY=ABS(V(ITJ)  J 
IF (VTRV.GT.VMAX) VMAXaVTRY 

110  CONTINUE  

CC  •• 

CC  ••  mRAPOLATE  U AND  V TO  EXTEBIUR  CINE 
CC  ••  .......ALONG  TOP 

00  112  liJiTK 

U(ltIYPU«3.*(U(I tIYPl-l»-U(I.lYPl-2) )»U(ItIYPl-3» 

V ( 1 i I YPi  j i3.*  rvTrtTYPi"=Tr=vTTrrYi»r=27TiTn;  lYPiO) 

CC  **  ....AND  ALONG  BOTTOM 

~ uTn  . * I iri  1 1 2 J I rt'3 ) I tj  Cl  ;tt 

v(iti)»3.*ivat2)-v(i.3n*vtifA» 

112  CONTINUE 

CC  *•  extrapolate  U ANO  V to  EXTERIOR  LINE 
CC  •*  .... . . . .'.  .lIL0fi6“WlSMT  'TTCE 

00  114  JsltlYPl 

irrTxpi.j»«3T»(uiixpi-rtj»-urixpr-2rjniU(rxpi-3tJ> 

V(IXPl.J)»3.*(V(IXPl-l.J»-V«IXPl-2tJ))»VtIXPl-3fJ) 

114  CONTINUE 

CC  ** 

IF  (KONT.LE.EOI  go  to ‘850 

CC  •*  TO  AVOID  SKIPPING 

■ 'TFIKONT.GEViri  T50  TO  '850 

CC  *• 

CC  •*  USE  EXTRAPOLATTOFTFOR  U ANO  V OW~THE~BOUNOARr 

CC  **  extrapolation  assumes  the  following  sequence  of  events 

CC  **  I SAVE  UIBliVIB) tAFTER-TmANGE-OF“REFERENCE  PT;)“-  “ 

CC  **  2 SAVE  U<B).V(B)  (AFTER  CHANGE  OF  REFERENCE  PT.) 

CC  •*  4 EXTRAPOCATE~'OECTA-CAPU 

CC  ••  3 EXTRAPOLATE  U(8»tV(8»  IF  DELTA  TIME  IS  SMALL  ENOUGH 

CC  **  5 IN  PLACC-OF^0REEN'*»-^WT-AtrTt«-rtBI— BY  EXTRnPOLATCO-OCcrr-CAPO 

CC  ** 

IF  (MOO(IGRNKTt3>-jE(>jOT-  GO-TO-BBO 

CC  *• 

K=M0O(  I0RNKTt31 

N0IF=1Y-2 

DO  760  I»2tIX - - - 


best  AVAllABlt  W 
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DO  760  J*2.IY.NOIF 
usn.jnrriirtTTjT 

VS(I«JtK)3V(ItJ) 

760  CONTINUE 
CC  ••  SAVE  time 

!F(k;CT.  n T5=TrRE 

CC  ••  SAVE  OELCU 

lF(K.£(J.I)-OEtCUS=OELCU 

CC  •• 

00  765  J32,IY 

US<lX»J»K)aU  (IXtJ) 

VS(IXtJfK)3V  lIXtJ) 

765  CONTINUE  

MRITE(  6<6026)  K 

6026  FORMAT  I"  U CB)  »V  <B)-  SAVED", IS) - 

CC  *• 

CC  **  

IF(M00(IGRNKT,3).NE.2)  00  TO  650 

CC  **  CACCUCATC  DELTA  TIME-  FOR  EXTRAPOLATION  USE  - - - 

nTIX»UMAX*0YI»VMAX»0XI*0T2 

OTCRIT»l./OTIX  --  - 

0T».8»0TCRIT 

TIMETaTlMe*-ai 

CC  **  ENO  CALCULATE 

N0lF=IV-2 

DO  780  l32«IX 

00  780  J>2«IX,N0IF 

C WRITE(6,6019)  I,US(IfJ,l) ,05(1,0,2) 

U ( I ,U)-»EXTRAP(US«1»0,1>  ,US-(  I-> 0, 3 ) , TS, T IME ,+mEX) 

VI 1,J)3EXTRAP(VS(I,J,1) ,VS(t,Jf2) ,TS,TIME  ,TIMET) 

C WRITE  (6,60l9)-iALC-l..J)-^V-LT,J) 

6019  FORMAT(I5,AF10.2) 

780  CONTINUE 

00  785  I»3,IVMI 

U( tX,  1 )jEXTBAP  IUS  1 1 X , I ,1 ) ,USHXr4 ,8)  ,T5,TfM&rT-»4€T-) 

V(IX,I)»exTRAP(VS<IX,I,I) ,VS(IX»I,2) ,TS»TIME,TIMET) 

C WRITE  16, 60191— I-,U+32-,-I4-^V-(-32vH 

785  CONTINUE 

WRITE (6, 6020) 

6020  FORMAT!"  UIB),V(B)  EXTRAPOLATED") 

DEL6U«6XTRAPI0EL-CU5,DEfcCUrT-5-r-T-IT4E ,TIHET) 

850  CONTINUE 

CC  

210  iprt*iprt*i 

C IFI IIPRT/IPRINT)»IP9INX,NE^IPRT+60-40-220 

CC  **  KONTF  * FINAL  COUNT 

CC  — •• - 

CC  *•  STORE  ON  TAPE  HERE 

C IF  IMOOIKONT,20).E»,0)-WRITE-HO)— STORE 

CC  ** 

CC  *•  CALL  PRINT  OUT-HERg- 

A :A)iAPRIlUTA 

CC  ••  - 

IFIKONT.EO. KONTF)  STOP  33 

220  CONTINUE — 

CC  •* 

IF(K0NT.6E.0)  GO-TO-100 ■ - - 

C IF (TIME.LT.TSTOP)GO  to  ioo 

CC  ••  — 

IFI  tIPRT/IPRINT)*IPRINT.EO.IPRT)GO  TO  225 

C WRITE  OUT  FINAL  ARRAYS- . 

CALL  PRTOUT 

GO  TO  (222,221,222,221)  ,ITYPE._  . 

221  CONTINUE 


BiSl  AVAILABLE  COPY 
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WHIIt(IF.lLC^liIORE 

2210  t<RITE(IPRTER«100«) 

1004  FORMATUH  t"-EIWAL-AaRAri_Q(sLQlSK__i_rEBMINATlON- NORMAL") 
STOP 

222  CONTINUE- - 

2220  NR1TE(IPRTER.1009) 

1009  FORMAT  tLM-U!_ 

STOP 

225  CONTINUE 


GO  TO  (2220t2210<2220t2210l tITYPE 
300  WRITE lIPRTERtlOOSJ 


1005  FORMATdH  t"TOP  LEFT  VALUE  OF  PSI  ON  BOUNDARY  MAS  CHANGED  TOO  MUCH 
1 EXECUTION  TERMINATED") 

CALL  PRTOUT 

- - STOP 

310  WRITE(IPRTER>1006) 

1006  FORMAT-d 

1 MUCH  EXECUTION  TERMINATED") 


CALL-  PRTOUT- 


STOP 

3000  CONTINUE 

WRITE dPRTERtlOOT) lERtARG 
1007  FORMATIIH  t*  ERROR  RETURN  FROM 
1 "ARO  . "fE16.8) 


r«T- 


END 

— subroutine  RRTOUT 


COMMON/BOX/OTtTIMEtOX>OY«HtPSI (33«63) • VORT (33t63f2) tU I33t63) t V (33t 

-l63^^^MA>TVMAX^eAPU^DELU^»1-^K0NT  

COMMON/BOXl/  FCFX(33t63) .FCFY«33.63) 

-C YVCITt33)B3) tVC2Tt33ia3) 

COMMON/lNDEX/lXtlY 

-COMMOlr^FCHEClt/  ICWTtIKNTttTCR 

IKMT»IKNT*I 

IXPI»IX»t 

IYP1»IY*1 


: PRINT  PSI-tUtVtVORT 

MRlTEIStlOOO) 

1000  FORMAT dMlT"PR-rirT- 


WRITE (6t999>TIME 
999  FORMAT  dH  -f»TTME- 
ISKIP-I 

M>l 


♦TCtGrBT- 


00  100  JsM.lYPltlSKIP 
L»IYP1— J»I- 


1001 


WRITE(6il001)L 

FORMAT  dHfr»"ROW  »-"t13T 

WRITE(6<1002) (PSI (ItL) tl^MtlxPlf ISKIP) 

-i»RITC(6tI0>e)  ) UfltL»Tl»H,|XPltlgHlP> 

WRITE  t6t  I 002)  ( Vd.L)«IaM.IXPltISKIP) 


HRITE(6tI002)  (FCFXdtL) 
WR I TE  ( 6t  1 0 02  MFGf^V-f  I-rtl-r 
WRITE(6il002)  (VClTdtL) 


tlaMtlXPltlSKIP) 
i- 


I tlaMtlXPltlSKIP) 
WRlTE«»rlOGa)  IVCef-d tL>  tlaN.l API, (SHIP) 


1002  FORMATdH  tI0E13.5) 
100  CONTINUE- 


IFdTER.GE.35)  STOP  2 
IF  ( IKNTjOEvKNTT-STGP-I- 
RETURN 

CNO 


BKf  AVAILABLE  COPY 


/kKED.SYStlB^  0l>-D5W5Y&l.TG0rLin,0tSt»-5MB 
//  00  OSN»FORT.SSPLI0«OISP«SHB 

//GOtSYSlfY  00 •- 


DEL  X«l/60  T|l)«TtU».Ol  CXPU«0»«1000 
-%  I NO  


tPRINT«20. 
ITAPE»V»  - 
TST0P«.1» 
PSCMT«.5» 


PSlTSTa.OOl tEPS«.001*RA«,25« 

— i6N^T«a»  ■ — 

ICNT*5. 


FITPX«l.»Fi^TPY»lT-i 

FITPX«1  .E50»FnPY»l  .E50t 



CAPUslOOO.t  . , , 

RUNN0«1 X ,9tFACT*.l «FACTX».303fFACT2»2.3A« 


K0NTF»A. 

(.ASTS>Ot 

IFlRST-l* 

1FIRST«0. 

SENO 


